FLOW ROUTING MODEL FOR UPPER INDUS RTVER (PAKISTAN) By DANIAL HASHMI B.Sc (Civil Engineering) University of Engineering & Technology Lahore. (Pakistan), 1984 A THESIS SUBMTTTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTERS OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1993 0DanialHashmi, 1993 In presenting this thesis in partial fulfillment of the requirements of an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Civil Engineering The University of British Columbia 2324 Main Mall Vancouver, Canada V6T 1W5 Date: ABSTRACT For the flow forecasting of the Upper Indus Basin (UIB) Pakistan, it is observed that large quantities of ungauged lateral flows enter many sections of the Indus River system during the monsoon and tiie snow melt season. The river is steep, flow travel times are short, and the data are currently available only on a mean daily flow basis. There are very few gauging stations and the travel time between them is less than a day, so that standard routing procedures such as the Muskingum method was difficult to apply because of the large amount of ungauged lateral flows entering the main channel. Calibrating the Muskingum model was therefore difficult because of the ungauged lateral flows, and also because the coefficients used in the Muskingum method are a function of routing period and vary considerably when the routing period is less or greater than the travel time. There was a need for a model which could not only route the flows tiirough the Upper Indus reaches but also could link all the sub basins contributing to the Indus flows. The whole UIB is divided into reaches, and a digital computer Link Model is developed which links these reaches and routes the flows through these reaches using the flow routing method developed by Quick and Pipes. The Link Model adds the known lateral flows to the main channel flows and also calculates the ungauged lateral flows if the stream gauging data for the downstream station is known. This Model can easily be interfaced with the UBC Watershed Model, which forecasts the land phase of the flow and the total system of models can be used for flow forecasting of UIB. u The routing method adopted is very flexible and easy to calibrate, because the principal feature of the routing method is that the routing coefficients are determined directly from the stage-discharge and area-discharge measurements at the gauging stations, so the method is pre-calibrated before commencing routing, and therefore avoids the problem of accounting for ungauged lateral inflows. The model has been tested on a 400 km segment of the main stem Upper Indus River, which lies above Tarbela dam i.e, between Besham Qila and Kachura gauging stations Fig (5.1). and was found satisfactory. Given the required data the Model can be used in any similar system. m Table of Contents Abstract ii List of Figures viii Acknowledgment xi 1 Introduction 1 2 Background and Literature Review 3 2.1 Introduction 3 2.2 Channel Routing Models 6 2.2.1 Empirical Models 6 2.2.1.1 Lag Model 6 2.2.1.2 Gauge relations 7 2.2.2 Linearized Models 7 2.2.2.1 Classical Wave Model 8 2.2.3 Hydrologic Models 8 2.2.3.1 Muskingimi Model 9 2.2.3.2 Muskingum Cunge Model 13 2.2.3.3 Kalinin-Miljukov Model 13 2.2.3.4 Lag and K Model 14 iv 2.2.4 Simplified Hydraulic Models 14 2.2.4.1 Kinematic Model 15 2.2.4.2 Diffusion Model 16 2.2.5 Complete Hydraulic Models 17 2.2.6 UBC Flow Model 18 Method of Analysis 21 3.1 Introduction 21 3.2 Theoretical Background 22 3.3 Model Description 29 3.3.1 Assumptions and Constraints for use of the Model 30 3.3.2 Channel Storage 31 3.3.3 Flow Translation 35 3.4 Flow Travel Time Test 46 3.5 Setup Procedure for the Flow Model 48 3.6 Flow Model Structure 54 3.7 How the Flow Model Works 59 3.7.1 Channel Routing 59 3.7.2 Lake Routing 61 3.8 Calibration of Flow Model 62 Link Model 69 4.2 Introduction 69 4.2 Link Model Structure 70 4.3 Setup Procedure for Link Model 73 4.4 How the Link Model Works 74 4.5 Link Model User Interface 76 4.5.1 Link Model 77 4.5.2 Basin Files 77 4.5.3 Graphics 77 4.5.4 Statistics 78 4.5.4.1 The Coefficient of Efficiency, E 78 4.5.4.2 The Coefficient of Determination, D 79 4.5.4.3 Volume Check 80 Study Area 81 5.1 Introduction 81 5.2 Hydrology of the Study Area 83 5.2.1 Indus River 83 5.2.2 Upper Indus Basin 83 5.2.3 Climate of the Karakoram 84 5.2.3.1 Precipitation 84 5.2.3.2 Temperature 85 5.3 Melt Process in Upper Indus Basin 86 5.3.1 Snow and Ice Redistribution 87 5.3.2 Glaciers as Stream Regulators 88 5.4 Flow Process in Sub-basins 88 5.5 Temporal Variability 90 5.6 Glacier Outburst 90 Application of the Flow Model to the Indus System 91 6.1 Introduction 91 6.2 Muskingum Model 92 6.3 Input Data 97 6.4 Discharge Pattern of River Indus 105 6.5 Link Model Results 108 6.6 Conclusions 111 References 122 VI Appendices 125 Appendix 1. Soiirce code for Link Model 126 Appendix 2. Description of Basin.dat file 142 Vll List of Figures 2.1 Modification in flood hydrograph brought by downstream movement of flood wave. 4 2.2 Prism and wedge storage in a channel reach. 9 3.1 Channel storage diagram for travel time TT = 1. 34 3.2 Typical water survey area-discharge curve. 36 3.3 A typical storage discharge curve and graphical representation of DSDQ values estimation for lake routing. 37 3.4 Kinematic translation of flood wave ignoring storage. 38 3.5 Effect of travel time on translated flows of River Indus at Partab bridge (1979 data). 41 3.6 Kinematic translation of flood wave ignoring storage for TT > 1 but TT< 2. 42 3.7 Kinematictranslationof flood wave ignoring storage for TT> 2. 44 3.8 Translation of inflow hydrograph due to low TT / DT values 47 3.9 Instabilities in routing due to large TT/DT values, low discharge values and sharp rise in inflows. 49 3.10 Regional map of the study area. 50 3.11 Typical water survey velocity discharge curve. 53 3.12 Constant data for three reaches main stem Upper Indus River. 56 3.12a Area-Discharge curve and DADQ values for River Indus at Kachura 57 Vlll 3.12b Velocity-Discharge curve and (DVDQ,QP,QQ) values for River Indus at Kachura 58 3.13 Recorded and calculated flow hydrographs of River Indus at Partab bridge. 63 3.14 Translation effect, this is controlled by changing DADQ values during calibration. 65 3.15 Storage effect, this is controlled by changing DADQ values during calibration. 66 3.16 Comparison of ungauged lateral flows between Besham and Kachura with a neighboring gauging station Astore. 68 4.1 Basin.dat file with job description prepared for Indus River between gauging stations Kachura and Besham Qila. 71 4.2 Main stem upper Indus River, schematic. 75 5.1 Map of the study area. 82 5.2 Flow chart for distribution ofwater input in sub basin. 89 6.1 Storage loop for River Indus reach between Kachura and Partab bridge using Muskingum method. 93 6.2 Recorded and Muskingum routed flows between Kachura and upstream of Partab bridge (1979 data). 95 6.3 Variationof Muskingum coefficients with K/DT values. 96 6.4 Schematic diagram of the study area showing available flow data duration. 98 6.5 Area discharge curve for River Indus at Kachura(l 991). 99 6.6 Velocity discharge curve for River Indus at Kachura(1991) 100 6.7 Areadischarge curve for River Indus at Shatial( 1991). 101 6.8 Velocity discharge curve for River Indus at Shatial (1991). 102 ix 6.9 Area discharge curve for River Indus at Partab bridge (1991). 103 6.10 Velocity discharge curve for River Indus at Partab bridge (1991). 104 6.11 Shift of rating curve for River Indus at Partab bridge. 106 6.12 Annual discharge hydrograph of River Indus at Partab bridge (1979 data). 107 6.12a Comparison of routed and recorded hydrographs of River Indus at Partab bridge (1979 data). 113 6.13 Comparison of routed and recorded hydrographs of River Indus at Partab bridge (1980 data). 114 6.14 Comparison of routed and recorded hydrographs of River Indus at Partab bridge (1983 data). 115 6.15 Comparison of routed and recorded hydrographs of River Indus at Partab bridge (1984 data). 116 6.16 Comparison of routed and recorded hydrographs of River Indus at Partab bridge (1985 data). 117 6.17 Comparison of routed and recorded hydrographs of River Indus at Shaital (1984 data). 118 6.18 Comparison of routed and recorded hydrographs of River Indus at Shaital (1985 data). 119 6.19 Comparison of routed and recorded hydrographs of River Indus at Besham Qila (1979 data). 120 6.20 Linear Plot of observed and calculated river Indus flows at Partab bridge (1979 data). 121 Acknowledgment I greatly appreciate the support and advice of my research supervisor, Dr. Michael Quick throughout my research and preparing this thesis. I also thank the other committee member Dr. Bill Caselton for his valuable comments on my thesis. Special thanks to my wife Sabiha and son Khurram for their continued support, and encouragement while this thesis was in preparation. I appreciate the help from Mr.Edmond Yu in improving my computer skills. This research is jointly funded by the Water and Power Development Authority (Pakistan) and the Intemational Development for Research Centre(Canada) through the British Columbia Hydro Intemational Limited in connection with the development of flow forecasting network for the Upper Indus Basin (Pakistan). XI Chapter 1 INTRODUCTION The objective of the study is to develop a flow routing method for a complex basin network so that flows can be forecasted at downstream stations with reasonable accuracy. The propagation of flow in space and time through a river or network of rivers is a complex problem. Flow propagation in natural rivers is complicated by several factors such as tributary inflows, variations in cross sections, variations in resistance as a function both of flow depth and location along the river, meandering of the river, interaction of the main channel and the flood plain and other various factors. A wide range of analytical methods is available, some quite empirical and others more physically based. Chapter 2 gives an overview of the various types of one-dimensional channel routing models. In this study the flow data used is from a mountainous area of the Upper Indus Basin (Pakistan). There are few gauging stations, the flow data available from these stations is on a mean daily basis and the flood wave travel time between these stations is usually less than a day. A conventional hydrologic method such as Muskingum method was not easy to apply due to the large amount of lateral flows entering the main channel and posing problems in calibrating the Muskingum model, and also because the Muskingum coefficients used in the model are a function of the travel time and vary considerably (Fig. 6.3) when travel time changes from a period less than to greater than the data interval. From the preliminary studies it became apparent that a nonlinear model was needed which could allow for the change of travel time with river stage and also with changed channel storage characteristics at different river levels. For this purpose the UBC Flow Model was tested and small adjustments were made to suit the extreme conditions. The results from this model were found satisfactory. This model is described in Chapter 3. There are many ungauged and gauged streams which enter the main channel at many locations. They constitute considerable volumes of flow and are active in different part of the year. To take care of these lateral flows and to route flows through the whole basin, a digital computer link model is developed. This Link Model includes a hydrograph combining component that performs the frmction of linking or combining the hydrographs from the various watersheds and tributaries. The frill description of this Link Model is given in Chapter 4. The physiography of the Upper Indus Basin is made up of a series of mountain ranges of extreme ruggedness and exceptionally high elevations, and these high mountains have a large impact on the weather system and on the stream flows. The flows in the study area are mainly due to the seasonal snow melt and glacier melt which are fed from the two principal sources of precipitation, most effective and dominant being the disturbances from west and the monsoons from the south. The weather effecting the study area is described in Chapter 5. Flow data from the study area was routed using the Muskingum method and the nonlinear UBC flow model. Results and performance of these models are presented in Chapter 6. Chapter 2 BACKGROUND AND LITERATURE REVIEW 2.1 Introduction :-Channel routing is a mathematical method (model) used to predict the changing magnitude, speed and shape of a flood wave as it propagates through waterways such as canals, rivers, lakes and reservoirs. The flood wave can emanate from precipitation runoff (rainfall or snow melt) or from reservoir releases (spillway flows or dam-failure). A flood hydrograph is modified in two ways as the storm water flows downstream. Firstly, the time of the peak rate of flow occurs later at downstream point. This is known as translation . Secondly, the magnitude of the peak rate of flow is dinainished at downstream points, the shape of hydrograph flattens out, and the volume of flood water takes longer to pass a lower section. This modification to the hydrograph is called attenuation and is illustrated in Figure 2.1. Routing methods are designed to estimate the modifications of the flood hydrograph as it moves downstream, so that flows can be predicted at downstream points. The theoretical foundation for channel routing was essentially achieved when Saint-Venant (1871) first derived the one dimensional equation of unsteady flow. The original Saint-Venant equations consist of conservation of mass equation: d(AV) ^ dA ^ ^ dx dt (2.1) and the conservation of momentum equation: dV .^^dV dh dt dx ^[dx ^ = 0 (2.2) in which t is time, x is distance along the longitudinal axis of waterway, A is cross sectional area, V is velocity, g is gravity acceleration constant, h is water surface elevation above datum and .S", is the friction slope which may be evaluated using a steady flow empirical formula such as Chezy or Manning equations. Upstream j r ' ^ ' " ' " ' " ° ' ' A Attenuated peak pownstream B FIG 2.1 MODIFICATION IN FLOOD HYDROGRAPH BROUGHT BY DOWNSTREAM MOVEMENT OF FLOOD WAVE. Due to the complexities of the Saint-Venant equations , their solution was not feasible, and various simplified approximations of flood wave propagation continued to develop. Simplified methods may be categorized as (i) purely empirical, (ii) linearization of the Saint-Venant equations, (iii) hydrologic; i.e. based on the conservation of mass and an approximate relationship between flow and storage, (iv) hydraulic; i.e. based on the conservation of mass and a simplified form of the conservation of momentum equation, (v) computer solutions of complete Saint-Venant equations. The choice of method depends very much on the nature of the problem and the data available. Flood routing computations are more easily carried out for a single reach of a river that has no tributaries joining it between the two ends of reach. According to the length of the reach and the magnitude of the flows considered, it may be necessary to assess contribution to the river fi-om lateral inflow, i.e. seepage or overland flow draining from, and distributed along, the banks. Further, the river network can be very complex system, and in routing flows down a main channel, the calculations are done for separate reaches with additional hydrographs being introduced for major tributaries, in order to develop flow forecasting for a major river system, detailed knowledge of main stream and various feeder channels is necessary. In this chapter a selection of flood routing methods will be presented, beginning with the simplest using the minimum of information and calculations progressing through to the more complex methods requiring the digital computer for their application. In the last classification and brief introduction of the channel routing model selected for the study area, i.e. the Upper Indus Basin, is presented. 2.2 Channel Routing Models 2,2.1 Empirical Models Some routing models are purely empirical, based on intuition and observations of past flood waves. Empirical models are limited to applications with sufficient observations of inflows and outflows of a reach of waterway to calibrate the essential empirical relationships or routing coefficients. They give the best results when applied to slowly fluctuating rivers with negligible lateral inflows and backwater effects. They require minimal computational resources; however, considerable effort can be required to derive the empirical parameters. Lag models and Gauge relations are examples of such models and are described below. 2.2.1.1 Lag Models Lag is defined as the difference in time between inflow and outflow within a routing reach. The successive average-lag method developed by Tatum (1940), assumes that there is a point downstream where the flow (/J at time {tj is equal to an average flow, i.e. (/,+/J/2. Tatum found that the number of successive averages occurring within a reach was approximately the time of travel of the wave divided by the reach length. Outflow (O) at the end of the reach is computed by: C*.. =C,/.+C,/,+... C^,/„„ (2.3) 7 where n is the number of subreaches (successive averages) within the routing reach. The routing coefficients used in the method can be obtained by trial and error using observed inflow and outflow hydrographs. The routing coefficients may also be obtained by a least squares correlation of inflow and outflow hydrographs, described by Linsley et al. (1949). A similar lag model known as the progressive average-lag method was reported by Harris(l 970). 2.2.1.2 Gauge relations These are pure empirical relationships which relate the flow at a downstream point to that at an upstream station, and are described by Linsley et al (1949). Gauge relations can be based on flows, water elevations, or a combination of each. The effect of lateral flows is automatically contained in the empirical relation. Gauge relation or stage routing usually relates the upstream and downstream stages by means of a correlation chart. A predicted flow at a downstream station can be readily obtained from the chart when the upstream stage or stages are given. 2.2.2 Linearized Models For theoretical and practical reasons, the full hydraulic equations of open channel flow are rarely used for channel routing, they are usually simplified to obtain solutions. The simplifications are either to totally ignore the least important nonlinear terms or linearize the remaining nonlinear terms in the equations. Given a sufficiently simpUfied form of the equations, they can be integrated analytically to obtain solutions of velocity and water surface elevations for any pair of {x,t) values with a relatively small computational effort. Usually the most common simplifying assumptions are (i) ignore the 8 second term in Saint-Venant momentum equation; (ii) constant cross sectional area, usually rectangular; (iii) constant channel bottom slope, often assumed to be zero; (iv) the friction slope term is linearized with respect to velocity and depth; (v) no lateral inflow; and (vi) the routed flood wave has a simple shape that is amenable to an analytical expression. An example is the classical wave model given below. 2.2.2.1 Classical Wave Model Neglecting lateral inflows, fiictional resistance and nonlinear terms V dAldx and V dVIdx in Equations 2.1 and 2.2 respectively, the following classical Hnear wave equations may be obtained: ^-SDj-r and ^ ^ ^ (2.4) where D is the average depth. The analytical solution to the equations (2.4) have the following form (Abbot, 1966): V =C:ix -4iDt) +q(x +Vi^) (2.5) where C/ and C/ are functions determined by initial flow conditions and boundary conditions. 2.2.3 Hydrologic Models This type of model uses the principle of continuity and a relationship between discharge and the temporary storage of excess volume of water during the flood period. The calculations are relatively simple and reasonably accurate and, from the hydrologist's point of view , generally give satisfactory results. Examples of this type of model are given in sections 2.2.3.1 to 2.2.3.4. 2.2.3.1 Muskingum Model In river routing the Muskingum method which was first developed by McCarthy(1938), is a commonly used hydrologic routing method for handling a variable discharge-storage relationship. This method models the storage volume of flooding in a river channel by a combination of wedge and prism storage as shown in Figure 2.2 below. r I -Q Wedge storage = / :X( / - (2) Prism storage FIG. 2.2 PRISM AND WEDGE STORAGE IN A CHANNEL REACH. 10 During the advance of a flood wave, inflow exceeds outflow, producing the wedge storage. During the recession, outflow exceeds inflow, resulting in a negative wedge. In addition there is a prism of storage which is formed by a volume of constant cross section along the length of prismatic channel. Assuming that the cross-sectional area of the flood flow is directly proportional the discharge at the section, the volume of the prism storage is equal to KQ where K is a. proportionality coefficient, and the volume of the wedge storage is equal to KX (1-0, where X is a weighting factor having therange 0 < r :^.5. The total storage is therefore the sum of the two components, S^KQ +KX(I-Q) (2.6) which can be rearranged to give the storage function for the Muskingum method. S=K[XI+{I-X)Q] (2.7) and represents a linear model for routing flow in streams. The value of X depends on the shape of the modeled wedge storage, and ranges from zero for reservoir- type storage to 0.5 for full wedge. In natural streams, X is between 0 and 0.3 with a mean value near 0.2. Greater accuracy in determining X is not done because the results of the method are relatively insensitive to the value of this parameter. The parameter K is flie time of travel of the flood wave through the channel reach. For hydrologic routing, the values of K and X are usually specified and assumed to be constant throughout the range of flow. The values of storage at a time J and j+1 can be written, respectively, as 11 SJ=K[XIJ+{I-X)QJ] (2.8) and Sj,, =^K[X I J,, -Kl -X) Qj,,] (2.9) using Eqs. (2.8) and (2.9), the change in storage over the time interval A ^ is Sj^i-Sj=K{[XIj*i+{l-X)Qf*i] -[XIJ+{I-X)Q]} (2.10) The change in storage can also be expressed, using Eq. (2.19). as S . , - s = ( i ± ^ A r - < ® ± e ^ A , (2,11) 2 2 Combining Equations 2.10 and 2.11 and simplifying gives fi/+i = Cilj^i + Czlj + CiQ (2.12) which is the routing equation for Muskingum method where Ci= ^f ^ f^ (2.13) C.=_4L^f^ (2.14) 2K{l-X)+At 2K{l-X)-At 2i5:(l-X) + A/ 12 where C\ + Ci +C^ =1 If observed inflow and outflow hydrographs are available for river reach, the values of A^ and X can be determined. Assuming various values of X and using known values of inflows and outflows, successive values of the numerator and denominator of the following expression for K can be computed. K .O-^^^U^'+^Mg'^'-^^)] n 16) X{lj^i+Ij)+{l-X){(3j^y-Qj) The computed values of the numerator and denominator are plotted for each time interval, with the numerator on the vertical axis. This usually produces a graph in the form of a loop. The value of X that produces a loop closest to a single line is taken to be the correct value for the reach, and K, according to Eq.(2.16) is equal to the slope of the line. Since K is the time required for the incremental flood wave to traverse the reach, its value can also be estimated as the observed time of travel of peak flow through the reach. The Muskingum equation has a constraint related to the relationship between the parameter K and the time interval At. Ideally the two should be equal, but At should not be less than 2KX to avoid negative coefficients and instabilities in the routing . 2KX <At <K (2.17) 13 To avoid this, a long routing reach is divided into subreaches so that the travel time through each subreach is approximately equal to the routing travel time. 2.2.3.2 Muskingum-Cunge Model A variant of the above method, sometimes referred to as the Muskingum-Cunge method, is detailed by Cunge (1969). The major feature of the Cunge modification is a means of estimating the parameters K and X and the time interval At from observed river geometry and flood wave behavior. 2.2J.3 Kalinin-MUjukov Model Another variation of Muskingum model is the Kalinin-Miljukov(1958) model. This model has the following form: O, = O, +(/. -O,) K, +(/, - / .) K, (2.18) where K.^l-e-'^"^" (2.19) K,=\-K^£ixl{c^t) (2.20) Ax=OISo{^yl AO) (2.21) in which A.ylAO is the slope of the depth-outflow rating curve. Equation (2.18) is identical to the Muskingum model if in the latter, K=Ax/c and X =0. 14 2.2.3.4 Lag and K Model Another hydrological storage routing model is the Lag and K model (Linsley et al, 1949). First, the inflow is lagged and then the outflow (O2) at time (^ 2) is determined by substituting the relation S2 -S\ =K (O2 -Oi) in equation: / .+ / , _ Q.+O, _ S,-S, At (2.22) and solving for 02,i.e.: 0^ = [l,+I,-0,{l-2K/At)]/{l+2K/At) (2.23) where K=Axlc (2.24) and c = ^ (2.25) dA The extent of lag and .^parameter may be constant, or they can be functions of inflow and outflow, respectively. 2.2.4 Simplifled Hydraulic Models The advent of computers made possible the development of another category of routing models based on Equation 2.1 and various simplification of Equation 2.2. These models are classified as kinematic, diffusion or quasi-steady depending on the terms retained in the Equation 2.2. Their applicability is mainly determined by the bed slope and wave shape. In general, the steeper 15 slopes associated with overland flow or steep streams with slow rising floods are amenable to use the kinematic models. The diffusion models have a wide range of applicability and can accommodate milder bottom slopes. However, there are still some practical combinations of mild slopping channels and flood wave shapes that are not suitable for either diffusion or kinematic approximations and are usually treated with complete Saint-Venant equations. 2.2.4.1 Kinematic Model There are many forms of the kinematic model; however the basic assumption used in their derivation is that equation (2.2) can be simplified and expressed in the following form: S^-S^ = 0 (2.26) in which dhl dx =dy/ dx—S^. Equation (2.26) implies that the momentum of the unsteady flow is assumed to be the same as that of the steady uniform flow as described by the Chezy or Manning equation or a similar expression in which discharge is a single-valued function of depth, i.e. dA/dQ =dA/dy =l /c in which Q is discharge (flow) and c is the kinematic wave speed and given by equation (2.25) Also, since dA _ dAdQ dt dQ dt ^- '^ and Q=AV, equation (2.1) can be written as the classical kinematic equation, i.e. 16 f . c | 2 = 0 (2.28) ot dx which can be solved by explicit or imphcit finite difference methods, the later being more efficient in most river apphcations. The kinematic wave model is limited to apphcations where single-valued depth-discharge ratings exists, and where backwater effects are insignificant, since in kinematic models flow disturbances can only propagate in the downstream direction. These models are very popular in apphcations to overland flow routing of precipitation runoff. Kinematic wave models have been used in channel routing by Harley et al. (1970) in the MIT catchment model and in the Hydrocomp model (Linsley,1971). A routine in HEC 1 package also deals with the kinematic method. When the kinematic approximation is used, flood wave attenuation, which is in fact often insignificant in actual steep streams, is not properly simulated because no dispersion terms exist; generally the simulated flood wave tends to change shape and steepen in the downstream direction. For the rivers with the bed slope of 0.001 or more the kinematic approximation often gives a fairly good representation of observed behavior. 2.2.4.2 Diffusion Model Flood waves in less steep rivers do not behave in a purely kinematic manner; the flood peak attenuates as the flood propagates downstream, and the stage discharge curve exhibits a looped behavior for rising and falling limbs of the flood hydrograph. Henderson(1966) examined this behavior in some detail. He writes the unsteady one-dimensional open channel flow equations in the 17 dx dt ^ ^ ^H ^2,dv^ ^2.30) dx g dt This can be rewritten in terms of slope ' ' - ' - - ^ - - ' ^ - - T ; (2.31) dx g ox g dt where H is tiie total energy of flow, v the velocity, B the surface width of the channel and U the local inflow. Henderson shows that flood wave attenuation can be represented in the dy momentum equation by including the surface slope term -~ to give dx 5- = S,-^ (2.32) d X When solved with the continuity Equation 2.29, the simplified momentum Equation 2.30 yields an expression similar to the heat diffusion equation. This model is tiierefore referred as the diffusion analogy ; it has been examined and explained by Hayami (1951), Henderson (1966). Ligget and Cunge(1975). Henderson (1963) shows tiiat the method can represent attenuation correctly and that the maximum discharge occurs before the maximum stage. 2.2.5 Complete Hydraulic Models 18 After the advent of high-speed computers much effort has been expended in the development of the complete (dynamic wave ) models. These are the solutions of complete Saint-Venant(l 871) equations and are categorized according to the method of solution, i.e, direct or characteristic methods. In the direct metiiod finite difference approximations for the partial derivatives are substituted dkectiy into Equations 2.1 and 2.2, and solutions are obtained for incremental times(A/) and distances (Ax) along the waterway. In the method of characteristics, the partial differential Equations 2.1 and 2.2 are first transformed into an equivalent set of four ordinary differential equations which are then approximated witii finite differences to obtain solutions. Dynamic models can be classified further as explicit or impUcit, depending on the type of finite difference scheme used in the solution. Explicit schemes transform the differential equations in to a set of easily solved algebraic equations. However, implicit schemes transforms the differential equations into a set of algebraic equations which must be solved simultaneously; the set of simultaneous equations may be either linear or non linear, the later requiring an iterative solution procedure. This type of model requires very detailed information on channel geometry and resistance. Usually the simple approximate models are adequate and are also more suitable for the requirement of most channel routing application. 2.2.6 UBC Flow Model The flow model presented in this thesis was developed by Quick and Pipes(1975) and is a modification of Lag and Route type model which can be 19 The flow model presented in this thesis was developed by Quick and Pipes(1975) and is a modification of Lag and Route type model which can be classified as a Simplified Hydraulic Model. In this model, effects of both kinematic wave and some diffusion are blended together to suit the steep river conditions, for which this model was developed. In this model the extent of lag and the storage parameter K both are variable and are fiinction of discharge. Storage parameter K is calculated using hydrometric stage discharge data while lag parameter TT is calculated using velocity-discharge relations. Routing is carried in two separate procedures, in the first procedure the inflows are routed through a simple linear reservoir using K = L * dA / dQ *STK (2.33) and storage outflow is calculated using QO(J) = QO(J-l) + [ QI(J) - QO(J-l) ] / ( 1 + K ) (2.34) where L is the reach length ,dA / dQ is the slope of the area-discharge curve and STK is the constant which defines the shape of the storage. In the second procedure these flows are translated by calculating travel time TT using Equation 2.35 and hydrometric velocity discharge curve TT = RTK/v (2.35) this travel time TT is used to calculate translated flows by using 20 Q0T(J+1) = QO(J) + (1 - TT) [ QO(J) - QO(J-l) ] (2.36) where v is the velocity of the flood wave read from hydrometric velocity -discharge curve, RTK is the inverse of kinematic wave travel time for the reach and QO(J)'s are the flows out of the channel storage calculated above. 21 Chapter 3 METHODS OF ANALYSIS 3.1 Introduction Channel routing procedures are used to calculate water levels and discharges along tiie river system during floods when unsteady flows in the form of "flood waves" moves through the channel system. Such waves can be many kilometers long and may be caused by changes in watershed flows lasting firom a few hours to several or many days. Flood waves are modified in the course of downstream travel. If there is no local inflow to the channel, the peak flow normally reduces in the downstream direction, a phenomenon referred to as flood wave attenuation. Because of hydraulic and sometimes mobile bed effects, the maximum water level often lags slightly behind the maximum discharge. (Henderson, 1966). The objective of the channel routing computations is to estimate flood wave travel time, modification and attenuation by accounting for the effects of the wave speed, channel storage and channel friction. A wide range of analytical methods are available, some quite empirical and others more physically based. The choice of method depends on the river channel characteristics, the accuracy required, the available data, and the time and budget available. 22 The flow model presented in this chapter is designed for the situation often encounter in the mountainous areas i.e, (1) The river slopes are steep to moderately steep. (2) The data base is sparse in which flows might only be available at daily intervals and at widely spaced locations in the channel system. (3) Very little information on channel geometry is available, and (4) Large amounts of ungauged lateral flows may come into the main channel. Application of the conventional metiiod is limited by the following factors: (a) The Muskingum method is not workable because the rivers are steep and the wave travel time between the gauging stations is less than the time interval of the flow data, thus resulting in negative, variable coefficients and instabilities in routing. These problems will be discussed later in Chapter 6. (b) Large amounts of ungauged lateral flows create considerable problems in determining the routing coefficients for the usual routing methods. (c) Although full hydraulic models can accurately simulate stream flow, these models require extensive data, which in most situations is not available. 3.2 Theoretical Background The Saint-Venant, 1871 Equations, 3.1 for continuity and 3.2 for momentum, are the governing equations for one-dimensional unsteady flow in open channels. 23 | £ + 5 i i = Si (3.1) dx dt | f+v|^+^|^-g(^.-&) = 0 (3.2) dt dx dx In which " y " = water depth; " B " = channel width; " x " = distance along channel;" i" = lateral inflow per unit length of channel; " So " = bed slope and " Sf " = frictional slope. These equations have various simplified forms, each defining a one dimensional distributed routing model. The momentum Equation 3.2 consists of terms for the physical processes that govem the flow momentum. These terms from left to right are: the local acceleration term, which describes the change in momentum due to the change in velocity over the time, the convective acceleration term, which describes the change of momentum due to change in velocity along the channel, the pressure force term, proportional to the change in water depth along the channel, the gravity force term, proportional to the bed slope" So " , and the friction force term, proportional to the fiictional slope" S/ ", The local and convective terms represent the effect of inertial forces on the flow. When a water level or flow rate is changed at a particular point in a channel carrying a sub-critical flow, the effects of these changes propagate both upstream and downstream. These backwater effects can be incorporated into distributed routing methods through the local acceleration, convective acceleration, and pressure terms. Different distributed flow routing models are produced using the fiill continuity equation but in which some terms of the momentum equation are eliminated. The simplest distributed model is the kinematic wave model, which neglects the local acceleration, convective 24 acceleration, and pressure terms in the momentum equation, which is a severe simplification, and reduces the momentum equation to. So = S/, so that the friction and gravity forces balance each other, an assumption which is very close to the truth for natural floods in steep rivers, whose slopes are of the order of 10 ft per mile or more. Relative magnitudes of these terms are given by Henderson, 1966, who writes the Equation 3.2 in terms of slopes as dy V dv idv So- — = S/ (3.3) dx g dx g dt and gives some typical values for tiie different slope terms taken from an actual river in steep alluvial country, which are ft/mile s ^y ^"' dx' 1 26, - , 2 V dv gdx ' 1 1 8 4 ' l ^ v gdt 1 20 (3.4) These figures relate to a very fast-rising flood in which the flow increased from 10,000 cusecs to 150,000 cusecs and decreased again to 10,000 cusecs within .24 hr. Even in this case, where the acceleration terms are comparatively large, the last three slope terms are so small that only "So " need to be retained in the momentum Equation 3.2. However, when the bed slope is very flat the dy/dx termmay well be of the same order as "So "and the significant term So-By Idx is to be retained in the momentum equation, which is usually referred to as the diffusion wave model. The fiiction term " 5/ ", is usually written in terms of Chezy" C " or Manning" n ". 25 Sf = T A ^ (3-5) Qn 149AR Sf = I , . . V „ a/3 (3.5b) For the kinematic wave model , discounting all the terms in the momentum equation (3.2) except " So " and using the Chezy form for"5/ ", Eq (3.2) reduces to For many natural channels," R " is approximately equal to the depth," y ", and the area can be written approximately as ( B.y ) . From Eq.(3.6) discharge Q is Q = CB/'^So"" (3.7) Differentiating with respect to " x" and substituting in Eq.(3.1) for (dQ/dx) (3.8) ,2 Jdx dt B From basic partial differential theory dy _ ^^y_ ,dy dt dt dx dt (3.9) Therefore the left hand side of Equation 3.8 will reduce to dyldt if dx/dt ={3ll)CSo^^ y''\ But from the Che2y formula, CS<,^^ y^^ is the channel velocity," v ", so that 26 — = 1 C So^'^y''^ (3.10) dt 2 reduces to " 3/2 v ". Eqs (3.8), (3.9), and (3.10) give the simple kinematic wave result that dy di . - dx 3 / - . U N -^ = — if = _ v (3.11) dt B dt 1 ^ ^ For zero lateral inflow 5/ this result shows that the depth" y ", should be a constant along the characteristic, defined by dxl dt =(3/2)v so that the flood wave will move through the channel system at 1.5 times the channel velocity, which is equal to the flood wave celerity. Linsley quote a range of flood wave celerities lying between 1.3 and 1.7 for different type of cross-sections and different resistance equations. There is no mechanism available in kinematic wave theory to subside or disperse the flood wave as it progresses through the channel. Subsidence may be controlled by the resistance and acceleration terms in the dynamic equations of motion, or by a much sin^ler mechanism, the pondage or storage introduced by the irregularities in a natural channel. These effects are dominant when the channel slopes are gentle. For the gentle slope case Q = BCy\y ^^ ' (3.12 ) dx by differentiating the equation (3.12) with respect to " x " gives 27 dl 2 dxT : > o - — — 2 ^x 7 y So dy dx (3.13) for the gentle slope case the approximate result as given by Henderson, 1966 of Equation 3.12 is 3 3 ^ (^ dy c = —V = —CJy\So—zr-2 2 r I dx (3.14) Substituting this in Equation 3.13 dx -1^ Bey d^ y dx J (3.15) For a kinematic wave with a wave velocity of c, it follows that for any point on wave profile dt dt dx (3.16) Eliminating dQIdx between Equation 3.1 with no lateral flows and Equation 3.14 leads to the result ^y - ^y 4^^y - ir^y _ _ _ __ rC—— — A i _ , dt dt dx dx^ (3.17) where 28 Kx = . "^ . (3.18) 3 V dx Comparing Equations 3.16 and 3.17 .The former is the standard wave equation common to all branches of mechanics and applicable in the present context to kinematic waves; it indicates that " y "appears constant to an observer moving with velocity c. Equation 3.17 also is a standard form of wave equation, containing an extra " diffusion " term K\^ yldx^ , which has the effect of making " y " decrease in the view of an observer moving with velocity c. This diffusion term includes the effect of the second slope term in Equation 3.4, which modify the kinematic character of the flood wave and makes it subside. Another source of subsidence could be channel irregularities which can behave as a chain of small lakes. In fact any form of off channel storage, such as seepage to the surrounding country will have the same effect. Hayami suggested that these storage effects could be indicated by a second diffusion coefficient Ki added to the first one, the combined coefficient K —K\ +Ki replacing K\ in Eq.(3.17). There is no theoretical basis for this proposal, but when applied in practice it gives satisfactory results even when K and c are assumed constant, which can only be approximately true. The great advantage of this " diffusion analogy " method is that if" K " and " c " are assumed constant, use can be made of a well-known wave equation which has known explicit solutions. Using these explicit solutions, the movement of the entire wave profile, however irregular it may be, can be traced with confidence. The form of the solution of the Equation 3.17 is very much dependent on the form of the initial flood wave introduced at the 29 upstream end of the reach; if, for instance the initial wave is of sinusoidal form then the Equation 3.17 has solution of the form y = yoe-" cos,{bt-ca) (3.19) indicating a progressing cosine wave moving at the kinematic wave velocity with an exponentially decaying amplitude, which is the classical solution to the diffusion problem. Flood waves are not usually simple cosine waves, but more realistic waves which could be represented by Fourier series. However, for real routing problems finite difference solutions are usually employed. 3.3 Model Description The Flow model is based on a simplified hydraulic equation of unsteady flow in open channel and it describes both the translation and attenuation effects of flood waves. In natural flood waves both kinematic and dynamic wave motions are present. However, in many cases, when the channel slope dominates, the flood waves behave approximately as a kinematic wave. Lighthill and Whitham (1955) proved that the velocity of the main part of a natural flood wave approximates that of a kinematic wave. Based on these ideas and realizing that the Eq. (3.19) can be decomposed simply into two effects i.e, a pure translation of flood wave and a decay function(subsidence effect). The Flow model is designed so that it independently models flood wave translation and flood wave subsidence in two separate procedures. The first procedure takes care of the decay function and is achieved by routing the inflows through a simple linear storage, which represent the storage 30 characteristics of the reach. The pure translation effect is achieved by the second procedure which mainly deals with the kinematic wave movement, it calculates the travel time from the Water Sxu^ey velocity-discharge curve, and using this travel time translates the outflow from the storage, through the reach. It is observed that this travel time is much more significant than the reservoir delay time. The routing coefficients for storage and translation are determined from the area-storage curves at the gauging stations, which are assumed to apply to the whole reach. This almost independent definition of travel time and channel storage has made the model very flexible and it is easy to fit it to a real system. Also, based on the theoretical groimds discussed above, the physical behavior of flood waves in a real channel system can be almost completely described by a travel time and a decay term, or subsidence, as modeled by a reservoir storage. Since the travel time and linear reservoir storages are related to such known characteristics of the river as area-discharge and velocity-discharge relationships, the model can be pre-caUbrated from the river gauging data and usually only small adjustments need to be made to these estimated routing factors. 3.3.1 Assumptions and Constraints for use of the Model Like most of the routing models, this model is based on the following assumptions. (a). The discharge data are given at equal time interval or routing period, and within this period the increase or decrease of inflow or outflow is assumed to vary linearly. 31 (b). The channel is divided into a number of reaches, each reach has practically constant characteristics. The flow is routed successively from reach to reach. In general, the reach length should be such that the average travel time of flood wave (TT), should nearly be equal to the hydrometric data duration (DT). For reasonable ftmctioning of Flow Model TT/DT should lie between 0.5 and 1.5 , but in any case TT/DT should not exceed the range 0.33 to 3.0 otherwise some serious instabilities such as negative routed flows can occur in routing. (c). Presently the model is working on a daily basis, accepting mean daily inflow data and the output is also on daily basis, but the model can easily be adjusted for any other inflow interval. 3.3.2 Channel Storage It is assumed that storage, S, and outflow QO relationship can be linearized for certain ranges of flow. Usually two or three linear segments give a close approximation to the actual curve. For a linear reservoir S = K{QO) (3.20) Continuity gives QI-QO = ^ (3.21) at For a linear process, each day's events are separable and independent of previous and future flows provided that yesterday's outflow is used as a reference level for both inflow and outflow [Quick and Pipes, 1975], i.e. 32 A{QI) = QI{J+l)-QO{j) AiQO) = Q0{J+1)-Q0{J) (3.22) (3.23) where " J " indicates the time period. The Equation 3.20 can be written as A{QI)-A{Q0) = f (3.24) Differentiating Equation 3.20 with respect to time will give dt At (3.25) Combining Equations 3.24 and 3.25 A{QI)-A{QO) = K ^{QO) At Rearranging the terms (3.26) A{QI) = A{QO) 14 At (3.27) Solving for A{Q0) and writing Kl At as a non- dimensional term K T A fee) = ,^^A{QI) (3.28) 33 By considering the daily flow increment as step increments and calculating each daily increment in inflow relative to yesterday's outflow, the complete reservoir routing procedure is described by combining Equations 3.12 and 3.27 A(fiO) = -1—[QI{J+i)-Q0{J)] (3.29) and the next days outflow is given by QO{j+l) = Q0{J)-^6.{Q0) (3.30) The change in storage is given by combining Equations 3.25 and 3.28 A5 = jf^ [QI{J+1)-QO{j)] (3.31) The procedure to carry out channel storage routing is so developed that it produces only a change in wedge storage, because the routing period equals the true travel time of the flood wave. The change in prism storage is taken care off by the translation process (to be discussed later). The value of storage parameter i^Jfor the wedge storage, Figure 3.1 is assume to be about half the value that would be calculated for a reservoir situation and this correction is made by introducing a constant STK, which is given the value 0.5 . The values of STK, which describes the shape of reservoir may vary according to situation, for very fast releases from reservoirs or for very fast rivers its value could be less than 0.5. QT(J} ^ ' INITIAL STEADY FLOW LEVEL ^ \ ^ QO(J) /Q0T(J+1) --^ .^.^ ^ Q0(J+1) WEDGE STORAGE ""'—••~..r""~~-~^ DEPTH ROUTING REACH FOR TT = 1 FIG. 3.1 CHANNEL STORAGE DL\GRAM FOR TRAVEL TIME TT = 1 35 The values of j^Tare derived directly from water survey data. Curves for area-discharge are plotted, as shown in the Figure 3.2 and the gradient is used to compute KT from KT = L — ^ ^ (3.32) dQ At in which L is the channel length. The corresponding relationship for calculating reservoir routing KT value is given by KS = —— (3.33) dQ At ' ^ in which dSldQ is taken from storage-discharge curve, obtained for a particular reservoir or lake. Figure 3.3. 3.3.3 Flow Translation The translation of the upstream hydrograph to some downstream point 1 day later is described by Figure 3.4, which is a plot of flow versus time. It is assumed that the upstream hydrograph is a linearly increasing flow. The flow is required for some intermediate downstream point which has a travel time, TT, from the upstream point. If the travel time were exactly equal to the routing period and the data interval DT, i.e, TT/DT = 1, then the downstream flow today would equal yesterday's upstream flow for daily flow data. If travel time is less than 1 as shown in Figure 3.4 , the downstream flow today will be greater then yesterday's upstream flow, i.e., discharge Z instead of discharge X. 20,000 8,000 0 50,000 100,000 150,000 200,000 DISCHARGE cusecs 250,000 300,000 FIG. 3.2. TYPICAL WATER SURVEY AREA - DISCHARGE CURVE ON w I (0 •D _1 I-O 111 9 UJ i CO 12.000 10,000 8,000 6.000 4.000 2.000 STRAIGHT LINE APPROXIMATION 0 2,000 4,000 6.000 8,000 10,000 . 12,000 14,000 16,000 DISCHARGE cusecs FIG. 3.3 A TYPICAL STORAGE DISCHARGE CURVE AND GRAPHICAL REPRESENTATION OF DSDQ VALUES ESTIMATION FOR LAKE ROUTING. FLOW TT(TT<1) yestardat's flow at upstream station TIME FIG. 3.4 KINEMATIC TRANSLATION OF FLOOD WAVE IGNORING STORAGE 00 39 Defining the upstream flow increment for the previous day as A g , in which AQ/At is assumed to be constant for each time period, the translated downstream flow, QOT{J +I), will be yesterday's inflow Ql{j) plus XZ. Thus QOT{J+I) = Ql{j)+XZ (3.34) By similar triangles ABC and CDE ZZ = {l-TT/DT) AQ (3.35) Therefore QOT{j+l) = QI{J)+{I-TT/DT)AQ (3.36) Equation 3.36 is used to translate the flow using the values of TT/DT calculated jfrom flow-velocity data for the reach and corrected for kinematic wave travel time by a factor RTK = v/c, whose values for different resistance equations and different channel cross sections are given by Linsley and Kohler, and ranges from 0.59 to 0.78 with the usually assumed values of 0.67. TT is calculated as TT/DT = RTK - — (3.37) V At if the values of TT/DT is close to one the Equation 3.36 works weU, but if the value of TT/DT is significantly greater than or less than 1 then the Equation 3.36 is modified according to the value of TT/DT. This modification changes the interpolation for AQ , and is described below for different cases. 40 Figure 3.5 shows a typical case of the variation of TT/DT with discharge in a steep river. Case I : When 1 < TT/DT < 2 Assuming flow to be linearly increasing, the translated downstream hydrograph is defined by Figure 3.6. Since the value of TT/DT is greater than 1 translated flow will be at point D instead of B, and the translated flow at point "B " will be the yesterday's upstream inflow, QI(J) minus" BD ". That is QOT{j+l) = QI{J)-AQ (3.38) From the Figure 3.6 ED = BC = (TT/DT - 1) ,comparing ABFG and A BED BD ED BG FG (3.39) QI{J)-QliJ-1) 1 ^^-^^^ Substituting the value of AQ from Equation 3.40 into Equation 3.38 and rearranging will give QOT{j+l) = Ql{j)4l-{TT/DT))[Ql{j)-QliJ-l)] (3.41) Case I I : When TT/DT > 2 41 ia /11 (Noiivana viva MOU) / (awii i3Avyi) 58 0\ a: CO <-> CM o m HI ft O 00 1 E M § i o Ed d E S 8 S soesno sMoid aaivisNvyi FLOW FIG. 3.6 KINEMATIC TRANSLATION OF FLOOD WAVE IGNORING STORAGE, FOR 1< TT / DT < 2 to 43 Following the same procedure as Case I, when the value of TT/DT is greater than 2, the translated flow QOT(J+l) will be at point " E " instead of point " I " on the flow-time diagram and it will be given by the following equation Q0T{J+1) = QI{J-1)-AQ (3.42) FromFig(3.7) CD = AD-AC = HE = TT/DT-2 Comparing AIEH and AIFG IE HE IF GF AQ _ {TT/DT)-2 QI{J -1) -QI{J -2) (3.43) (3.44) Substituting the value of A g in the Equation 3.42 and rearranging will give Q0T{J+1) = QI{J-I)-^2-{TT/DT))[QI{J-I)-QI{J-2)] (3.45) Case III: When TT/DT < 0.50 This condition is mainly associated with high flows and it means that the flows are faster by a factor two or more than anticipated. (Ideally the value of TT should be equal to the time interval of the flow data At ,ie, TT/DT =1). Normally imder the condition when TT/DT < 1, Equation 3.36 is used to calculate translated flows. In this equation, the change in inflow, AQ, is TRAVEL TIME TT FLOW FIG. 3.7 KINEMATIC TRANSLATION OF FLOOD WAVE IGNORING STORAGE AND FOR TT / DT > 2 ^ 45 calculated by extrapolating the upstream flow data for fiiture days. In order to limit the extrapolation and for the stability of the flow model (which will be discussed later in section 3.4 ) the reach length is selected in such a way that TT/DT values should always remain above 0.33. To bring TT/DT values further close to one and to bring compatibility between (1 - TT/DT) and the extent of extrapolation for calculating AQ, time steps in upstream flow data interval are introduced. These time steps are introduced for two future days, because only two days upstream flows can effect downstream translated flows. Therefore, two upstream inflows QI(J-K).5) and QI(J+1.5) are calculated by assuming linear increase in upstream flows and by averaging QI(J); QI(J+1) and QI(J+1), QI(J+2) respectively. Routing equations for storage and translation are adjusted accordingly by revising the value of routing period from DT to half DT , then calculating KT by using KT = RL * DADQ * STK * 2 / DT (3.46) and then using Equatons, 3.29 and 3.30 to calculate outflows from the storage. These outflows are then translated by using following three equations: V = (L0G{Q0{J +15)) -(LOG{QQ)))*DVDQ (3.47) TT/DT = RTK*RL*2/{V*DT) (3.48) QOTiJ+2) = Q0{J+15)+{I-TT/DT)(Q0{J+2)-Q0{J+15)) (3.19) where QO(J)'s are the out flows from the storage corresponding to QI(J)'s. 46 Although the cases described above present the solutions to the particular practical routing problems, usually the reach lengths are selected in such a way that the bulk of the flows have a travel time close to one, and within the limit of 0.5 TT/DT to 1.5 TT/DT. Nevertheless there are flows which lie beyond these limits, and therefore the Flow Model has been tested for the extreme cases and is described below. 3.4 Flow Travel Time Test Movement of the flood wave has a significant effect on the shape of the wave profile. In the Flow Model, flood wave movement is measured by the wave travel time (TT) over a fixed distance i.e, reach length. It is desired to keep this travel time equal to the inflow data interval (DT), i.e, TT/DT=1. This TT/DT has a considerable bearing on the translating part of the Flow Model. For this reason, and to see how far the model can be stretched, the model was tested for a large range of different values of TT/DT, brought about by varying reach lengths. The findings from the tests are presented below: For TT/DT « 1 For very low values of TT/DT i.e (0.25 to 0.06) a hydrograph for upper Indus River was generated (see Fig. 3.8 ). Two phenomenon are evident from these hydrographs, first the translation effect can be seen in the complete shifting of tfie calculated hydrograph in the direction of flow. Secondly very Uttle steepening of the wave can be seen in one of the peaks. Figure 3.8 may suggest that ultimately a surge will form. However^even for quite low values of o (0 o LU O OH < o CO Q 140,000 120.000 100.000 80.000 60.000 40.000 -20.000 •v...../. TT/DT LEGEND INFLOWS TRANSLATED FLOWS >-/>^ .^--' J_ 0.25 0.2 0.15 0.1 0.05 0 30 60 90 120 150 180 210 240 270 300 330 360 TIME days Q FIG 3.8 TRANSLATION OF INFLOW HYDROGRAPH DUE TO LOW TT / DT VALUES 4^ 48 TT/DT i.e 0.05 to 0.25, no surge was observed. In order to prevent any chance of kinematic shock the lower limit for TT/DT value has been fixed as 0.33. For TT/DT » 1 For high values of TT/DT ranging from 1.5 to 10 the inflow hydrograph for the above river reach has been routed to the downstream point. For TT/DT values greater than 3.0, coupled with the steep rise in inflows, serious instabilities in the form of negative flows were developed in tiie routing, particularly in translating the flows through the river reach. This is shown in Figure 3.9 From the above results, it would seem that a good working range for TT/DT values is 0.33 to 3.0 and these values can be used when establishing a suitable reach length. 3.5 Setup Procedure for the Flow Model Before running the Flow Model, the following preliminary steps are completed 1. A Regional map is prepared which demarcates the basins and sub basins in the region and which clearly shows the position of gauging stations, the distances between them and points at which lateral flows enter the main channel. Figure 3.10 is a regional map from the Upper Indus Basin showing the above features. 2. Data collection: (a) Historical stream flow data is obtained, (b) Flow characteristic data such as, stage-discharge , area-discharge, and velocity-discharge values for a wide range of flows at each gauging station on the 200,000 150,000 o 100,000 UJ O o Q 50,000 (50,000) NEGATIVE FLOW 20 40 60 80 TIME days FIG. 3.9 INSTABILITIES IN ROUTING DUE TO LARGE TT / DT VALUES , LOW DISCHARGE VALUES AND SHARP RISE IN INFLOWS 100 UJ i I 5 i LU ai > 4^ i ' U O f A » ( , O FIG. 3.10 REGIONAL MAP OF THE STUDY AREA 51 channel is obtained. Table (3.1) is a typical sheet showing characteristics of flow data at a gauging station. 3. The channel is divided into a suitable number of reaches. The reach lengths are determined so that the bulk of the flows through the reach should have travel times close to the time interval of the flow data, and the beginning or end of some of the reaches should coincide with the gauging station. 4. DVDQ values: For each channel reach, the mean flow velocity and the corresponding discharge data from the upstream or downstream gauging station is obtained (sometimes the average of bodi gauging stations is used). Velocity is plotted against the log of discharge, and this velocity-discharge curve is divided into approximate straight line segments as shown in Figure(3.11). Slopes of each segment are determined, and these are the DVDQ values corresponding to each range of flows. The discharge value at the junction of two segments is the QP value for the lower segment. The discharge values obtained by extending segments to the abscissa are QQ values corresponding to that range of discharges lying under that segment. These values are required to calculate velocity as a function of discharge in the model. 5. DADQ values: The same procedure as described in (4) above, is adopted, but this time cross sectional channel flow area is plotted against the corresponding discharge values, these values are obtained from the representative gauging station located either upstream or downstream of the reach. Slopes of the segments obtained from the straight line approximation are the DADQ values and junction points of the segments are QA values. These values are used in channel storage routing calculations, where the storage parameter KT is calculated by using Eq.(3.32) 52 Discharge measurements of Indus River Near Besham Qila 1991 1 2 3 4 5 6 7 NO 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Date 21.1.91 11.2.91 20.3.91 9.4.91 29.4.91 8.5.91 10.5.91 27,5.91 27.5.91 27.6.91 08.7.91 13.7.91 22.7.91 05.9.91 12.9.91 width Ft 430 438 434 454 452 464 466 479 474 575 577 582 597 574 570 Area Sq.ft. 4450 5170 5310 6340 6410 7190 7530 10500 8920 17100 18000 21800 21900 15600 14500 Mean velocity Fps 3.71 4.7 5.18 5.52 6.24 9.21 6.95 11.43 11.04 12.4 12.5 12.98 15.8 11.73 11.36 Gauge height Feet 12.19 14.1 14.36 16.34 16.7 17.6 18.51 24.6 21.35 37.58 39.15 45.84 50.27 34.45 32.5 Discharge eft. 16500 24300 27500 35000 40000 66200 52300 120000 98500 212000 225000 283000 346000 183000 165000 rABLE 3 . 1 TYPICAL FLOW CHARACTERISTIC DATA AT A GAUGING STATION 16 14 12 ^ 10 8 > 0 Straight line approximations Pivot points 9.37 9.57 9.80 9.88 10.22 LN (DISCHARGE) cusecs 10.81 11.29 11.65 FIG. 3.11 TYPICAL WATER SURVEY VELOCITY - DISCHARGE CURVE 54 6. DSDQ values: In order to route the flows through a lake in the river channel system, a storage versus the lake outflow discharge curve is defined. In the moimtainous areas where generally the lakes are steep sided, sufficiently accurate storage can be estimated as S = A{Z-ZO) (3.50) where zo is the elevation of the bed of the outlet or stream, or other convenient datum, z is the lake surface elevation and A is the surface area of the lake. Flow values from the outflow gauging station of the lake is obtained and a plot is drawn between storage and discharge. The same procedure as described above in (5) is adopted to calculate DSDQ segment slopes, which are shown in the Figure(3.3). These values of DSDQ are used to calculate routing parameter KS from Eq.(3.33) 3.6 Flow Model Structure The flow model is basically designed for computations by computer and is coded in Quick Basic 4.5 . The main feature of the program is its modular structure, in which the program is broken into simple, manageable and conceptually independent tasks, called procedures. These procedures are easy to understand, maintain and debug, and they can be defined once and yet executed any number of times throughout a program. The model is controlled by a user friendly interface which describes all steps and data requirement for running the flow model. When the model has been run there is provision to view the output from the model in tabular or in graphic form, and user can compare graphically the routed flows with that of the recorded flows. 55 Statistical analysis of the measured and calculated hydrographs is also available, and this analysis produces measures of the model performance, expressed in terms of an efficiency based on residual variance, and voliraietric accuracy. Data requirement: Two types of files are required by the Flow Model (a) Inflow file : This type of file contains historic stream flow data at various gauging stations, such as upstream inflow data for routing and downstream flow data for calibration of the flow model. Each gauging station has a separate file and these files are named to briefly describe the flow measuring station with an extension of FLW, for instance, the data file for the gauging station on the River Indus at Partab bridge can be named as IND@PRB.FLW. To give fast computer access a binary file structure is used, witii each record containing the date of occurrence and the corresponding flow. Because of ease in handling, dates are expressed as Julian dates . These files can be prepared from the standard water survey data files by using a conversion program (MKBRYFLW.EXE). (b) Constant file: This ASCII file contains reach characteristic data. Usually a gauging station at the upstream end of the reach adequately describes the reach characteristics, but sometimes an average of upstream and downstream or downstream values are judged to be more representative. Data included is reach length, time interval of flow data, DVDQ, DADQ, DSDQ values etc. One set of data is required for each reach. Fig (3.12) shows a typical set of constant data information required to run the flow model. Figure 3.12a and 56 *** RIVER INDUS CONSTANT VALUES RL,L,M,DT,QO(0),RTK,STK BETWEEN KACHURA AND PARTAB BRIDGE *** 420000.00 4 QQ,QP,DVDQ 500.00 21000 6900.00 40000 8100.00115850 53100.00300000 QA, DADQ 22000.00 40000.00 116000.00 4 86400.00 400.00 0.670 0.500 00 00 00 00 1. 2, 2, 9. 410 076 500 375 300000.00 0.164 0.122 0.086 0.019 ***RIVER INDUS CONSTANT VALUES BETWEEN PARTAB BRIDGE AND SHAITAL RL,L,M,DT,QO(0),RTK,STK 400.00 0.670 0.500 5 86400.00 .00 ,00 00 00 0. 2. 7. 8, 753 818 727 409 492000.00 4 QQ,QP,DVDQ 500.00 27200 13400.00 80800 40150.00156400 47100.00400000 QA, DADQ 8200.00 0.683 48000.00 0.036 115000.00 0.043 160000.00 0.036 400000.00 0.025 *** RIVER INDUS CONSTANT VALUES BETWEEN SHATIAL AND BEESHAM QILA *** RL,L,M,DT,QO(0),RTK,STK 365000.00 2 QQ,QP,DVDQ 500.00 21600.00 14000.00118000.00 50000.00400000.00 QA, DADQ 22000.00 0.450 93000.00 0.069 170500.00 0.044 400000.00 0.017 4 86400.00 400.00 0.670 0.500 0, 2. 7. 570 764 250 FIG. 3.12 CONSTANTS DATA FOR THREE REACHES MAIN STEM UPPER RIVE INDUS zr in 2i Q 14.000 12,000 10,000 8,000 6,000 4,000 2.000 DADQ = 0.086. QA = 40000. -QA =2200C DADQ = 0.122 DADQ = 0.164 I I QA= 116000 DADQ = 0.0193 I 0 20,000 40,000 60,000 80,000 100.000 120,000 140,000 160,000 DISCHARGE cusecs FIG 3.12a AREA DISCHARGE CURVE AND DADQ VALUES FOR RIVER INDUS AT KACHURA 14 12 10 S. o DVDQ = 2.076 QQ2 = 6900 QQ3 = 8100 - / 115850 /DVDQ = 9.375 QQ 4 = 53100 8.5 9.5 10 10.5 Ln DISCHARGE cusecs 11 11.5 12 12.5 FIG 3.12b VELOCITY DISCHARGE CURVE AND (DVDQ, ,QQ, QP) VALUES FOR RIVER INDUS AT KACHURA oo 59 Figure 3.12b illustrates the method adopted to arrive at these constant values only for the first reach between Kachura and Partab bridge. The information for the preparation of the file can be obtained fi-om a regional map, area-discharge curves, velocity-discharge curves, storage-discharge curve and the inflow data. This file can be prepared with the program (CONSTMKF.EXE). 3.7 How the Flow Model Works 3.7.1 Channel Routing To run the routing program, the user must supply the basic information, such as inflow file name, constant file name, number of reaches and starting and finishing dates for running the model. After checking the acceptable format and presence of the input data, the program creates the working inflow file by extracting the specified duration of inflow data firom the total historic record inflow file. Then the program opens the constant file and reads the data required for storage routing. A loop is started and every flow is routed through the linear storage using equations(3.51) and(3.52). KT/DT= L * DADQ * STK / DT (3.51) and storage routed outflow is calculated as follows QO(J) = QO (J-1) + [ QI (J) - QO (J-1) ] / ( 1+ (KT/DT) ) (3.52) 60 In order to translate these outflows, QO(J), through the reach, another computational loop is started. First the velocity is calculated from the velocity discharge relationship by using v = [ Ln (QO(J) - Ln (QQ) ] * DVDQ (3.53) in which QQ is the x-intercept value of the segment drawn for straight line approximation on the velocity discharge curve, and DVDQ is the slope of that segment, read from the constant file. The ratio of the travel time to the time interval of inflows is calculated from TT/DT = RTK / ( v * DT ) (3.54) in which RTK is the inverse of kinematic wave speed constant and RL is the reach length. The value of TT/DT is used to translate the flows from the storage through the reach and is also used to determine the best set of inflows for estimating the downstream flows. The selection of the best equation for translation is carried out through a series of IF statements, as follows: Case 1 when TT/DT = 1 Q0T(J+1) = QO(J) (3.55) Case 2 when TT/DT >1 but TT/DT <2 Q0T(J+1) = QO(J) + (1-TT/DT) * [QO(J) - QO(J-l)] (3.56) 61 Case 3 when TT/DT >2 QOT(J+l) = QO(J-l) + (2-TT/DT) * [QO(J-l) - Q0(J-2)] (3.57) Case 4 when TT/DT < 1 but TT/DT > 0.5 Q0T(J+1) = QO(J) + (1-TT/DT) * [Q0(J+1) - QO(J)] (3.58) Case 5 when TT/DT < 0.5 A time step is introduced in the inflow data by taking an average of the inflows, these flows are routed from the linear storage using half value of DT for calculating KT/DT i.e. KT/DT = RL * DADQ * STK * 2 / DT (3.59) and translating these flows through the reach by using Equations 3.60 and 3.61 TT/DT = RTK * RL * 2 / (v * DT) (3.60) Q0T(J+1) = QO(J+0.5) + (1-TT/DT) * [Q0(J+1) - QO(J+0.5)] (3.61) where QO's are the storage outflows. 3.7.2 Lake Routing 62 If any lake is present in the system routing through is carried out in the same manner as the storage routing, except that the routing coefficient KS/DT is calculated from Equation 3.62, where DSDQ values are the slopes values read from the storage-discharge ciu^e, as shown in Figure 3.3 KS/DT = DSDQ / DT ( 3.62) and the routing is done by using Equation 3.52. There is no translation of the flow in lake routing because the travel time in a lake is normally very fast, and can be ignored on a daily basis. 3.8 Calibration of the Flow Model For the Flow Model, we can derive values of the routing coefficients by using the standard water survey data, which is measured to establish the stage-discharge relationship. The required data is the mean channel velocity and cross sectional area as a function of discharge. If this data is representative for the reach the model converges to the recorded downstream flows at the first attempt, and no calibration is needed. There are channels for which some refinement of the coefficients is needed and an iterative procedure is carried out as described below: The model is run and the behavior of the routed flows are seen by comparing them with the recorded flows at the downstream gauging station. This is done by plotting the routed and calculated hydrographs on the same scale and on the same graph as shown in Figure 3.13, Important features to note are (a) the general trends of both the hydrographs. (b) the shift (lags) and 250,000 200,000 8 150,000 0) lU O < <J 100.000 50.000 I I I I LEGEND RECORDED FLOWS CALCULATED FLOWS ± 0 30 60 90 120 150 180 210 240 270 300 330 360 TIME days FIG. 3.13 RECORDED AND CALCULATED FLOW HYDROGRAPH OF RIVER INDUS AT PARTAB BRIDGE ON 64 points of shifts in the hydrographs. (c) the peaks of hydrographs, and (d) the rising and falling limbs of hydrographs. As discussed above the changes which can occur in the shape of the hydrograph due to the routing are the subsidence of peaks, horizontal shift and the steepening of the wave front. Any other effect is regarded as the influence of the ungauged lateral flows. The horizontal shifts or steepening of wave front which is caused by the wave velocity and is the characteristic behavior of the kinematic wave is controlled by DVDQ values. If the calculated hydrograph is shifted in the flow direction, then it indicates that the flows are moving faster and tiiey can be controlled by reducing the DVDQ values . Figure 3.14 indicates this behavior. The storage component of the model introduces some dispersion and attenuation and retards the steepening of the wave front of the routed hydrograph. Figure 3.15 indicates the same behavior, in Has figure significant attenuation of peak flows and very little delay in steepening of the flood wave is evident. This type of situation is controlled by changing DADQ values. DADQ values inversely effect the routed flows, if the routed flows are significantly less than the recorded flows then DADQ values are decreased. After adjusting these values the model is run again and the routed and recorded flows hydrographs are compared. If any significant difference is observed the above procedure is repeated. Ungauged lateral flows create considerable problems in calibration. Usually these ungauged flows obscure any clear trends of translation error and attenuation error. When these ungauged lateral flows are considerable, the routing model is used to estimate them by subtracting routed upstream flows from the recorded downstream flows. To do this, the model coefficients must be assumed to be correctly determined from the water survey data, because no w o W O LU O < X o CO Q 140.000 120,000 100.000 80,000 60,000 40,000 20,000 0 LEGEND RECORDED FLOWS PRECAUBRATED TRANSLMED FLOWS ± ± 0 30 60 90 120 150 180 210 240 270 300 330 360 TIME days FIG. 3.14 TRANSLATION EFFECT, THIS IS CONTROLLED BY CHANGING DVDQ VALUES DURING CALIBRATION J5 o </) O in O < X o w Q 140,000 120,000 100,000 80.000 60,000 40,000 20,000 LEGEND RECORDED FLOWS PRECALIBRATED ROUTED FLOWS 30 60 90 120 150 180 210 240 270 300 330 360 TIME days 0 \ FIG. 3.15 STORAGE EFFECT, CONTROLLED BY CHANGING DADQ VALUES DURING CALIBRATION ^ 67 calibration is possible when there is large ungauged lateral flows. Then the estimated lateral flows are compared with other neighborhood recorded data, if any is available. For example, while routing flows through three reaches between Kachura and Besham Qila considerable volumes of ungauged lateral flows enter these reaches. These flows were calculated by subtracting the routed flows from the recorded flows, then these flows were compared with the nearby watershed Astore flows, lying in the same basin. In order to find the degree of statistical hydrologic homogeneity between the ungauged lateral flows and Astore flows correlation and linear regression tests were carried out. The results of both the tests showed that there is a strong correlation between the flows with the sample correlation value of 0.81 and the standard error of lateral flows estimate as 13419 which is about 10% of the maximum flows. Figure 3.16 illustrates the same comparison, showing that the overall trend, distribution of flows and timings of the storm events or melt response from both the basins are nearly similar, therefore the difference of recorded and estimated flows at Beesham Qila are truly ungauged lateral flows. 140,000 120,000 100,000 g 80,000 Ui o a: < o 60,000 w Q 40,000 20.000 30 ASTORE FLOWS' 25 UNGAUGED LATERAL FLOWS 60 90 120 300 330 150 180 210 240 270 TIME days FIG. 3.16 COMPARISON OF UNGAUGED LATERAL FLOWS BETWEEN BESHAM AND KACHURA WITH A NEIGHBOURING GAUGING STATION ASTORE (1980) 360 00 69 Chapter 4 LINK MODEL 4.1 Introduction: The Link Model is used to connect all the tributary watershed inputs, the channels, sub-reaches, lakes and reservoirs into a total main channel system. It is a menu driven computer model, coded in Quick Basic 4.5 in which all the principal features contributing in the channel routing are divided into components. One of the first steps in developing inputs to the Link Model is the preparation of a basin layout or schematic diagram that can be used as a guide in organizing the data. In this basin layout, all the important features of the watershed and channel system are indicated on a map of a suitable scale. These features are the hydrologic boimdaries of sub-basins, the main channel, position of lake (if any), diversion and confluence points of tributaries and gauging station locations on the tributaries and the main channel. In the preparation of the routing scheme the main channel is divided in to a number of reaches, end points or meeting points of channel are node points, these node points are given numbers to identify different reaches in the network. This is easily carried out by schematic which helps clarify the correct sequence of sub-basin runoff, routing, and the flow combining computations. The routing is carried out by a routing component which is fully described in Chapter 3. The hydrograph combining component, essential to the 70 operation of the system as a whole, performs the function of linking or combining the hydrographs from the various watersheds and the tributaries. This linking process is repeated until the whole river system or sector of river system is included into the network. This chapter describes the basic structure of the Link Model, the components of model, their function and the data involved to operate them, the steps involved in preparation of input data, and the different facilities avaUable to view, print and compare the routed and recorded flows at gauging stations. The source code for the Link Model is given in Appendix . 1. 4.2 Link Model structure The Link Model is a computer program which is written in structured Basic and is divided into six procedures. Each procedure is conceptually an independent task and describes a physical phenomenon occurring in the basin. Prior to the running of the model the whole basin is divided into a nimiber of sub-basins and the river channels into a number of reaches. In order to carry out different tasks in a specific order, a set of instructions is passed to the main model through BASIN.DAT file. The program follows the sequence of the instructions, which invoke the appropriate procedures. These instructions themselves follow the sequence of physical phenomenon occurring in the basin. Figure 4.1 shows a typical set of instructions for the Link Model and following is the description of each procedure: 1. NEWREACH : This is usually the first instruction to the program and it describes the start of a new reach. Whenever this information is given to the program it adds one to the previous reach number and these reach numbers are 71 INSTRUCTIONS AND INDUS 4 INDSCNST.DAT IND@KCH.FLW NEWREACH CHANNEL ADDLATL GIL@ALM.FLW FILEOUT NEWREACH ADDLATL AST@DyON.FLW CHANNEL FILEOUT NEWREACH CHANNEL FILEOUT JOBEND INFORMATIONS FOR BASIN MODEL FOR UPPER INDUS BASIN Ncune of the channel Number of nodes File name for constant data Inflow file name Start of first reach Perform channel routing Add lateral flows at this point Lateral flow file name Prepare output files at this node Start of second reach Add lateral flows Lateral flow file name Perform channel routing Prepare output files at this node Start of third reach Perform channel routing Prepare output files at this node Finish the program FIG. 4.1 BASIN.DAT FILE WITH JOB DESCRIPTION PREPARED FOR INDUS RIVER BETWEEN GAUGING STATIONS KACHURA AND BESHAM QILA. 72 used as an extension to the working files. The results generated from each reach is stored in these working files, and for later use, the output from any reach is identified by the reach numbers. 2. CHANNEL ROUTING : This procedure performs the channel routing and is completely described in Chapter 3. 3. LAKE ROUTING : This procedure performs lake routing by treating it as simple linear level pool reservoir, and is described at section 3.7.2. 4. ADD LATERAL: This procedure combines two hydrographs at common outlet point of sub-basin and channel. It opens the latest channel flow file and the lateral flow file, reads the date of the first entry in the channel flow file and then traces the matching date for the lateral flow, after getting the dates matched it simply adds the corresponding flows and store them in the latest channel flow file. 5. FILE OUT : This procedure prepares output files at the downstream end of the reach. Files are usually named after the name of the channel and the extension to the files is the node point i.e., the point where usually the output is required. These files are further required to view the results, draw hydrographs, compare the results or generate statistics. 6. JOB END : This is the last instruction to the model and it terminates the routing process and switches to the main menu. 73 4.3 Setup Procedure For Link Model Before running the Link Model, the necessary input data is prepared. This data may include the inflows in the beginning of the channel reach, lateral flows either from the gauging stations or outflows from the sub-basins watersheds and their point of confluence with the main channel, and the reach characteristic data. All this necessary data is first collected and then organized in different files as described in the following steps: 1. Preparation of inflow files: The inflows to the channel system must all be specified before routing can be carried out. These tributary inflows will correspond to the various gauging stations throughout the system. In addition to this there will probably be some inflows that are not measured, and these will have to be estimated. For model calibration, the gauged inflows will be available as measured data. The ungauged flows may have to be estimated by routing the known upstream flows and then subtracting these routed flows from the downstream recorded flows. For forecasting the flows in the main channel, the various tributary inflows will be estimated with the watershed model and they will be available as lateral flows to the main channel. In order to keep the data condensed and make the program efficient, the flow files are kept as binary. The file structure for both type of flows i.e, main channel, tributary or lateral wiU remain the same. Each record is of eight bytes, with first four bytes for date(in Julian calendar) of occurrence of flow and rest for the corresponding flows. 74 2. Preparation of the constant file: These are ASCII files and they carry information on the reach characteristics of the channel. Preparation of these files are explained in the Oiapter 3, with an addition that information on each reach is arranged in the sequence as it appears in the real river system. If there is any lake in the system then it will be treated as a single reach and the constant values such as DSDQ, QS etc., will appear on the constant file accordingly. 3. Preparation of the Basin.dat file: This file is in ASCII mode and carries the (i) Information which describes the basin features, (ii) The series of instructions which control the routing and linking of flows and (iii) the name of different files, tiiat are used at different stages of program. The First few lines of the file describes the file description, channel name, node numbers, constant file name and inflow file name, and rest of the file describes the series of instructions to the main program. These instructions are passed on to the main program which performs the linking and routing operations and is described below. Basin.dat file can be prepared from BASINMKF.EXE program. Figure. 4.1 is an output of the Basin.dat file prepared for the main stem of the upper Indus River and Figure. 4.2 illustrates schematically the main stem of upper Indus River. 4.4 How the Link Model Works The Link Model follows series of instructions described by the Basin.dat file, all these instructions activate certain tasks which are defined above in Model structure. Basin.dat file also includes the names of different LEGEND GAUGINGSTATION RIVERS RESERVOIR SHAITAL - © t o : ALAM 0 BRIDGE PARTAB BRIDGE^ KACHURA RIVER INDUS FORT BEESHAM DVON UJ <tt: TARBBELA DAM FIG. 4.2 MAIN STEM UPPER INDUS RIVER, SCHEMATIC 76 files which are used at different stages of routing, such as the file containing the characteristic data for different reaches, and the files containing tributary inflows at different node points. Figure 4.1 is a sample of the Basin.dat file prepared for the upper stem of tiie River Indus, and in this file all the instructions for routing and the linking of flows are prepared firom the schematic diagram shown in Figure 4.2. Prior to the running of the main program the user selects the Basin.dat file and then enters the period of run. The main program opens the selected Basin.dat file and reads the instructions in sequence and perform the task accordingly. Program also produces different files at the node points which are used later on for different types of outputs, this is explained in the next section. Detail description of the instructions laid down in the Basin.dat file, which was prepared for the Upper Indus River is given in Appendix 2. The independent definition of different tasks has made the Link Model very simple, easy to understand, interpret and easy to fit to tiie real situation. The key to the Link Model operation is the preparation of Basin.dat file which is explained in Appendix 2, and Figure 4.2. 4.5 Link Model User Interface The User interface is tiie term applied to the parts of the software that enable the user to interact with the program. These include the screens that appears while the program is running and the keys used to move around the screens, enter data, select from lists and menus and perform similar operations. There are input screens that enable the user to give instructions and information to the computer, and output screens that display information derived by running the program. 77 In the Link Model, the input screens have a title at the top to indicate which function it is a part of, and the help line at the bottom to provide the information about the keys and the messages from the system. The main menu screen which comes up at the start of program has a bar menu at the top that lists the available Link Model operations and, in the help line at the bottom of the screen, a description of the menu option at the current cursor position. Following is the description of the items available on the main menu 4.5.1 Link Model: The main program that runs tiie Link Model for the selected basin and selected period of time. 4.5.2 Basin Files: Displays the names of the files describing the information of the basin, names of the files carrying inflow and lateral flow data and the operating sequence instructions to the main program. 4.5.3 Graphics: Many of the characteristics of the estimated and recorded flows can be viewed through the graphical presentations. Following are some of the plots available for comparing and calibrating the Flow model: 1. Calculated flows 2. Inflows to the reach, flows routed from the storage, and tianslated flows. 3. Recorded and estimated/calculated flows. 78 4. Error signal, this helps in examining the distribution of and spotting extreme values of estimated residuals. 5. Accumulated errors, this is important to evaluate the long term effects of estimated residuals. 4.5.4 Statistics: The overall performance of an estimated hydrograph can be evaluated against the observed hydrograph using a statistical approach which is based on the sum of the squares of residuals. The following criteria are available to assess the performance. 4.5.4.1 The Coefficient of Efficiency, E The coefficient of efficiency , E , can be expressed as follows (Nash and Sutcliffe, 1970) E = ^ (4.1) where S = ^{Qo{t) -Qe{t)y (4.2) " • -^M So = ^{Qo{t) -Qo) (4.3) where Qo is the observed flow; Qo is the mean of the observed flow for the given record Qe is the estimated flow. 79 The quantity S is the direct measurement of the discrepancy between the observed hydrograph and the estimated hydrograph. A high value of S indicates a large discrepancy and therefore a poor routing, while a small value indicate the opposite. So is the value of S if the calculated hydrograph is a horizontal line equal to the mean recorded flow. The coefficient of efficiency, E, is the difference between So and S normalized to So. Therefore, E is a non-dimensional measure of the quality of the estimated flows. A value of E close to unity indicates a good estimation, while a bad flow estimation is characterized by an E value much less than unity. 4.5.4.2 The Coefflcient of Determination, D The coefficient of determination, D, can be expressed as follows: D = ^ (4.4) where Sr = I too -Qr{i)f (4.5) where Qr is the flow calculated from linear regression of the recorded flow against the estimated , so that the criterion D measures the performance of the estimated flows afl;er removing residuals linearly related to the estimated flows(usually called linear errors). It is also an indicator to determine the correlation between the observed and calculated hydrographs. A very small coefficient of determination D, means that the two hydrographs are only weakly correlated, thus indicating error in the data, missing information on laterals or a serious problem in model structure . 80 4.5.4.3 Volume Check: A volume check is an important criteria to assess the model performance, one of the basic equations used in the development of Flow Model was the conservation of mass Equation 3.1. A final test of this equation can be seen by calculating the difference of the total volume of the water to the channel and the total volume of water from the channel. Ideally this should be zero but due to retention of some volume in storage, assumptions made in extrapolating and interpolating flows and round-off values during calculations, there may be some small volume error. Provided this error is small, mass is essentially conserved and this is important especially when estimating ungauged lateral flows. 81 Chapter 5 STUDY AREA 5.1 Introduction: All the test data used in the Flow Model and Link Model are taken from the area drained by the Upper Indus River between the gauging stations Kachura and Beesham Qila and is shown in Figure 5.1. The study area is located in northern Pakistan and encompasses an area of some 160,000 sq. km in a moxmtainous region of the Western Himalaya and Karakoram Ranges, and is a basin having an extreme degree of ruggedness, range of relief and numerous glaciers. There are two principle soiirces of precipitation, westerly disturbances which usually occur in winters and monsoons from south east in late July and August. Climate and topography have a great impact on the magnitude and timing of flows in the river Indus. River flows are composed primarily of water from snow and ice melt, while the monsoons may be important in the more southerly basin. The importance of these three sources varies through time. This chapter deals with the general description of the area, the major inputs contributing to the flows in the River Indus and the governing mechanism behind them. FIG. 5.1 MAP OF THE STUDY AREA oo 83 5.2 Hydrology of the Study Area 5.2.1 Indus River The Indus River has its headwaters in the high arid Tibetan Plateau area of China, and flows westward through Ladakh Province of northem India into northem Pakistan, (see Fig. 5.1). The river continues to flow to the west past Skardu, around Nanga Parbat (El. 8125m), and then flows south to Tarbela and on to the Arabian Sea. The total drainage area of the Indus River just upstream of Tarbela is 162,000 sq. km, and the average annual flow at Beesham Qila is 2350 ciunecs. The main tributaries of the Indus River in the study area are the Shyok, Shigar, Hunza, Gilgit and Astore. 5.2.2 Upper Indus Basin The area and the other physical processes in the Upper Indus Basin are described by Hewitt (1988) and are presented below: The Upper Indus Basin consist of a series of moimtain ranges of extreme ruggedness and exceptionally high elevations. The ranges run approximately west-east in an arc from the Hindu Kush and the Pamirs in the west and north to the Himalayan chain to the south and the east. The latter culminates in the Nanga Parbat massif, which forms a barrier to the northward movement of monsoon storms. The Central Karakoram forms the major part and water sources of the Upper Indus. With hundreds of peaks in excess of 6000 m elevation, it consist of a wide belt of mountain ridges and high valleys. On the greater Karakoram range there is extensive formation of glaciers. This is in contrast to some of the lesser ranges, and much of the greater Himalayan chain to the east whose high, but already not very extensive, massifs are 84 separated by deeply incised antecedent rivers, and do not support glaciers of such extent or concentration as those in the Karakoram. The glaciers of Karakoram cover an area of about 16000 sq. km and it is the ablation zones of these large glaciers which produce the bulk of glacier ice melt. 5.2.3 Climate of the Karakoram 5.2.3.1 Precipitation There are two principal sources of precipitation in the Indus Basin. The most important source of precipitation is from low pressure systems coming from west and this is the major moisture source during the winter months. A secondary summer accumulation season at high elevation in the Karakoram and Hindu Kush is usually derived from the same source. These westerly sources provide the dominant nourishment for the glacier system of the Karakoram and provides very important flows in the months of April, May and June and contribute significantly in July and August. The second source in terms of total amounts of precipitation delivered is the monsoon, and these storm systems come from the south. These storms are restricted in time from late July through September and are generally limited to the Main Himalayan and Front Ranges. Mayewski and Jeschke(1979) and Mayewski et al.(1980) have shown that in some years the monsoon is strong enough to break through the Front Ranges and can deliver substantial precipitation to the central Karakoram region. This was also confirmed by Wake (1987), who showed that there were substantial snowfalls in the central Karakoram derived from the summer monsoon. 85 A common characteristic of the study area is great spatial variability in precipitation over quite short distances and also regionally. The Karakorams exhibit exceptional spatial and especially altitudinal variability . Two strong regional trends may be seen. As a result of the monsoon influence, the Front ranges towards the south and the east of the basin receive heavy precipitation. Farther north and west this source of moisture is less significant. The westerly depressions give most of their precipitation to the western and northern portions of the region while there are noticeable gradients towards the center and east of the basin reflecting shadowing by the major ranges. The combination of these trends produces extremely high precipitation on the main Karakoram Range, semi-desert conditions occur in much of the interior valleys, plateau and ranges of the eastem Upper Indus and Shyok watersheds. Extreme local relief is associated with extreme altitudinal variability in precipitation. The deep and often extensive valley bottom are very dry, while the high mountains close by receive heavy precipitation. The bulk of aU precipitation occurs in an altitudinal zone roughly between 2500 m and 6000 m elevation. The elevation of maximum precipitation as given by Hewitt (1988) is probably about 5000m in the westem Karakoram, 6000m in the eastem Karakoram. 5.2.3.2 Temperatures Surface air temperatures show extreme spatial variability and depend primarily on elevation, but with strong seasonal and diumal cycles. Solar radiation , the prime source of energy for the melting of snow and ice over most of the basin, also shows great variability as a function of elevation, surface slope, aspect and changeable weather conditions. North facing and 86 south facing slopes commonly have quite different radiation regimes, producing marked difference in local hydrology. The high temperatures and very dry atmosphere normally prevailing throughout most of the summer within the valley system of the Upper Indus Basin promote high evaporation rates. In many of the low lying areas within the Karakoram evaporation is so high that surface runoff only occurs during rare, torrential thimder storms. As a result of the combination of low precipitation coupled with intense evaporation in these valleys bottom situation, water yield from these areas below 2500m is negligible. 5.3 Melt Process in Upper Indus Basin In describing the hydrology for any mountain range it is convenient to think of rainfall, snow melt and icemelt as the three major components for water yield. Each of them is space and time dependent. All are dependent on elevation, season and location within the region. In the Karakoram freezing conditions prevail for part of the winter at aU elevations, so that nearly all the precipitation falls as snow. Thus most of the precipitation derived from westerly depressions is snowfall. In July and August, temperatures are above freezing up to 5000 m. Hence most of the monsoon precipitation falls as rain, especially in the Front Ranges. While rainfall makes its way quickly into the river system (with a lag dependent on surface and subsurface conditions), snowfall remains in place (apart from the local redistribution by winds, avalanches and glacier flow) until the spring and summer. The snow melt yield from lower elevations is modest and rapid, occurring in early spring . Yield from progressively higher elevation is late in the summers, compressed into fewer months of the year, but is greater in total because of the much deeper 87 snowpacks. Glaciers represent long term storage reservoirs of water and glacier melt from the Karakoram region and other glaciated basins is considerable during the melt season. 5.3.1 Snow and Ice redistribution Hewitt(1988) describes that, there are very significant transfers of snow and ice within the hydrological system prior to melt. Snow is redistributed after falling by two major processes. Wind action redistributes snow locally (usually within a few hundreds of meters of where it originally fell). Deflation in some areas accompanied by drifting and great thickening in adjacent areas occurring during winter has important implications for the details of snow-melt processes taking place in the spring and summer seasons. Avalanching of both snow and ice is a major process of redistribution in the Karakoram. It greatly effects the nourishment process of most glaciers in the region and has major implications for melt hydrology in that vast quantity of snow and ice are transported to much lower elevations very quickly (primarily during spring time but continuing at high elevation throughout the summer) where they melt out relatively slowly but consistently. Huge quantities of ice are steadily transported downslope as a result of glacier movement from the accumulation areas, where snowfall exceeds melt to the ablation areas where there is net annual melting of ice. Compared to avalanching or melt water flow, glacier movement is slow. It may take hundreds of years for individual particle of ice to move from place of original deposition to the place of melt, but it is a constant process. Karakoram glaciers, unlike most continental glaciers, are relatively vigorous, 88 and the latest studies by Hewitt (1985) indicate flow rates in the middle ablation areas in excess of 150 m per annum. 5.3.2 Glaciers as stream regulators One of the main characteristics of glaciers in the Upper Indus Basin is to regulate stream flow. This characteristic as explained by Bradely (1990) is a consequence of radiative exchanges at the glacier surface associated with the snow cover. In cool wet summers ample precipitation falls as rain at low elevations increasing river flow, while at higher elevations snow fall is more likely which increases the surface albedo and therefore increases the reflectivity. It is estimated that for fresh snow fall this reflectivity is 90%. i.e only 10% of energy is available for melting. This loss of energy will therefore prevent melting. In hot summers, however a high input of radiation and absence of fresh snow increases seasonal melting of glacier ice which can supplement river flow. The importance of ice-melt is thus enhanced at such times when total flow from runoff of rainfall on valley sides is little. 5.4 Flow process in sub-basins The daily amounts of snow melt, icemelt and rain are subject to various delays and loss before eventually appearing as river flow at some watershed outflow point. The snow melt, icemelt and rain inputs are subdivided between evapotranspiration loss (which includes evaporation and sublimation loss), ground water zone, subsurface zone and impervious area according to the soil and weather conditions, and emerge as a fast, medium and slow flows . This whole process is described in Fig .(5.2.) RAINFALL.SNOWMELT & ICEMELT EVAPOTRANSPIRATION EXCESS OfSIRII^ WAte? 2DI*E IMPERVIOUS AREA -^^ FAST FLOW O — ^ * SLOW FLOW SUB SURFACE ZONE MEDIUM FLOW TOTAL FLOW FIG. 5.2 FLOW CHART FOR DISTRBUTION OF WATER INPUT IN SUB BASIN 00 90 5.5 Temporal Variability Temporal variability in discharge is another outstanding feature of Upper Indus Basin. Summer snow melt in combination with glacier ice melt produce the very high discharge peaks of mid and late summer. In the Front Ranges, however, spring and summer snow melt is later augmented by monsoon rains to produce even higher instantaneous discharges. 5.6 Glacier Outburst Glacier outburst floods occur throughout the Karakoram and several classic events have taken place this century; documentation of these has been given by Hewitt (1982). 91 Chapter 6 APPLICATION OF THE FLOW MODEL TO THE INDUS SYSTEM 6.1 Introduction The Basic purpose of the study is to develop some mechanism which can easily and effectively route the flows from the whole river system to some downstream point of the basin. To test the flow modeling, the Upper Indus Basin (Pakistan) was selected and the hydrometric data from the sub-basins and the main channel for the upper stem of River Indus between Kachura and Partab bridge was used. Analysis of the data showed that the stream flow records are scattered, so that only seven years of continuous data is available for the area under study. These inflows are on a daily basis and the flood wave travel time between gauging stations is less than a day, in addition to this there are significant amoimts of ungauged lateral flows entering the main channel. For comparative purposes, the flow model method was compared with the widely used conventional Muskingum method. It was found that the application of the Muskingum method is difficult, because the ungauged lateral flows entering the main channel posed problems in calibrating the Muskingum model. The flow routing method proposed by Quick and Pipes(1975) as described in Chapter 3, was found satisfactory for flow routing. In this model lateral flows are not a serious problem, because the model coefficients are pre-determined and can be calculated from water survey data, and allowance can 92 be made for the variation of travel time with discharge The Flow model can be used to determine these ungauged lateral flows when they are not directly measured. A Link Model has been developed for the study area in which the hydrographs from various watersheds and tributaries are combined and the routing is carried out by using the above flow model. This Link Model also calculates the ungauged lateral flows. In this chapter first the application and limitation of the Muskingum method is discussed then the problems encountered in the use of the input data are given. The effect of climatic changes on the River Indus is described. The results generated after running the Link model are discussed and finally the conclusion from the study is drawn. 6.2 Muskingum Method The Muskingum method was applied on upper Indus River reach between Kachura and Partab bridge (see Fig 4.2) for the flow data of 1979. For the inflows the data at the Kachura were used and for the outflow data the recorded flows at Partab bridge minus the recorded flows of River Gilgit at Alam bridge were used. In order to calculate the Muskingum routing coefficients first K and X are estimated, for this storage loops calculations were carried out using the method proposed by Chow(1964). It was found that some ungauged lateral flows were entering in the reach resulting in negative storage values. This resulted in difficulties in loop formation. Various lateral flows values were tried, for 8% of lateral flows a loop was formed, then different X values were tried to collapse the loop (see Fig. 6.1) but little progress was made, approximate K values were estimated from the loop and routing coefficients were calculated. Then these coefficients were used to route ID CO O H 00 140,000 120,000 100,000 80,000 60,000 40.000 20,000 K = 35/30 =1.16 days x = 0.2 20,000 40,000 60,000 80,000 X I + ( 1 - X) O cusecs 100,000 120,000 140,000 FIG. 6.1 STORAGE LOOP FOR THE RIVER INDUS REACH BETWEEN KACHURA AND PARTAB BRIDGE USING MUSKINGUM METHOD ^ (>i 94 the Kachura inflows to a downstream point just above River Gilgit and River Indus confluence. Routed flows were compared with recorded flows minus estimated ungauged lateral flows and are presented in Figure 6.2. It is evident from the Figure 6.1 that this scheme for estimating routing coefficient does not conform to the conventional method of forming a narrow loop. Flow travel time for the same data was calculated using kinematic wave estimates and it was found that there is a large variation in the kinematic wave speed. The K/DT values varied between 0.35 and 1.3. This was compared with the Muskingum coefficients values for different K/DT values calculated from the following : „ X-2x(KIDT) 2{KIDT){\-X)+1 ^ ^ ^^ ^ \+2x{KIDT) 2{KIDT){\-X)+1 C3 = 2 ( ^ / Z ) r ) ( l - A r ) - l 2{KIDT){l-x)+l ^ ' ^ where DT is the flow data duration and x is assumed as 0.2. These Muskingum coefficient values for different K/DT are plotted and is shown in the Figure 6.3. It is evident that for K/DT values less than one there is a variation in the coefficient values, which may results in instabilities in calculations. 14Q00O MUSKINGUM ROUTED FLOWS, i w CO p 1201000 1001000 aaooo eotooo 4Q00O 201000 150 165 180 I 196 210 TIME days 225 240 255 FIG 6.2 RECORDED AND MUSKINGUM ROUTED FLOWS BETWEEN KACHURA AND UPSTREAM OF PARTAB BRIDGE (1979 DATA) 270 ai I CO O CM o O UJ o u. u. UJ o O (D 2 5 O z « 3 8 (0.5) h K/DT VALUES FIG. 6.3 VARIATION OF MUSKINGUM COEFFICIENT WITH K / DT VALUES ON 97 6.3 Input Data Input data has an important role in setting up, calibrating, and evaluating the performance of the Model. Inadequacy of input data may undermine the quality of the flow forecast for one or more of the following reasons: Insufficient data record length : Model parameters calibrated from short data record are unreliable. On the other hand, long data records provides the opportunity to detect systematic errors in the data. Fig (6.4) gives the overall picture of the available years of record for the flow data used in the Flow Model. It is clear from Figure 6.4 that only two years of flow data at Shaital gauging station is available to run Link Model and check the system imder study as a whole. However, the Link Model can be run for seven year period i.e, 1979 to 1985 and the results can be compared at Partab bridge and Besham Qila. There are different sources of errors in the flow measurement, such as human or equipment errors, which can effect the reliability of the data. During the analysis of the area-discharge and velocity-discharge ctirves from the study area, for almost all the gauging stations some errors were observed These curves are used in calculating the channel routing parameter for storage and flow travel time for translation and are shown from Figure 6.5 to Figure 6.10. Common observations to these curves are that they are not rising gradually, there are some points on the curves which describes the abnormal behavior such as decrease in flow area with increase in flow discharge and decrease in average flow velocity with increase in discharge. These curves were corrected GILGIT AT GILGIT (1966 TO 1971) () INDUS AT SHATIAL (1984 TO 1987) HUNZAATDANYOR (1975 TO 1981) $ GILGIT AT ALAM BRIBGE (1979 TO 1985) kO ASTOREATDOYAN^ (1979 TO 1988) INDUS AT FARTAB BRIDGE (1979 TO 1986) RIVER INDUS . (^ INDUS AT KACHURA (1979 TO 1986) il w X w a. INDUS AT BEESHAM QILA (i979T01987) LEGEND GAUGING STATION RESERVOIR STUDY AREA BOUNDARY TARBELA DAM FIG 6.4 SCHEMATIC DIAGRAM OF THE STUDY AREA SHOWING AVIALABLE FLOW DATA DURATIONS 14,000 12.000 -20,000 40,000 60,000 80.000 100,000 120,000 140,000 160,000 DISCHARGE cusecs FIG. 6.5 AREA DISCHARGE CURVE FOR RIVER INDUS AT KACHURA (1991) 14 12 -10 -O Q 6 -4 -20,000 40,000 60,000 80,000 DISCHARGE cusecs 100,000 120,000 140,000 160,000 FIG. 6.6 VELOCITY DISCHARGE CURVE FOR RIVER INDUS AT KACHURA (1991) o o 22,000 20,000 18,000 -t 16,000 i 14,000 12,000 -10,000 -8,000 0 50,000 100,000 150,000 200,000 DISCHARGE cusecs 250,000 300.000 FIG. 6.7 AREA DISCHARGE CURVE FOR RIVER INDUS AT SHATIAL (1991) o 50,000 100,000 250,000 150,000 200,000 DISCHARGE cusecs FIG. 6.8 VELOCITY DISCHARGE CURVE FOR RIVER INDUS AT SHATIAL (1991) 300.000 o 16,000 14,000 -12,000 a-S 10,000 8,000 6,000 0 50,000 100.000 150,000 DISCHARGE cusecs 200,000 250,000 FIG. 6.9 AREA DISCHARGE CURVE FOR RIVER INDUS AT PARTAB BRIDGE (1991) o ^ ^ ^ 50,000 100,000 150,000 200.000 250,000 DISCHARGE cusecs FIG. 6.10 VELOCITY DISCHARGE CURVE FOR RIVER INDUS AT PARTAB BRIDGE (1991) o 105 before DADQ, DVDQ etc., values were determined from them. These corrections for Partab bridge data is shown in Figiu-e 3.12a and Figure 3.12b Equipment and human factors are not the only causes of flow measurement errors, since flow rating curves may not represent all of the flow conditions especially during high flood events due to scour of the channel section, or following the flood events due to the sedimentation in the channel at the gauging station. This phenomenon was observed for the flow data at Partab Bridge gauging station where the rating curve was shifted in one year and is shown in Figure 6.11. 6.4 Discharge Pattern of River Indus The annual recorded discharge hydrograph of River Indus at Partab bridge for the year 1979 is shown in Fig.(6.12). The form of the recorded hydrograph reflects the nature of radiative exchanges between snow and ice-melt and the effect of weather pattem on the flows. Discharge remains low over the winter, from January to February, during which time the majority of precipitation is stored within the snowpacks, only ground water and snow melt at lower elevations supplements river flow. Discharge first increases in mid March when there is a small rainfall-induced runoff peak. From March the flow of Indus rises sharply in a series of stages to the peak flow in June as snow melts from progressively higher elevations in the main ranges. Rise of the snowline exposes an increasing area of glacier ice, melting of which is normally responsible for subsequent peaks in discharge after June. The form of these peaks is effected primarily by the passage of reoccurring weather systems over the mountains which deposit fresh snow on 70 66 60 ti 55 S3 in i 50 w 1991 RATING CURVE 45 40 1992 RATING CURVE 35 501000 1001000 1501000 DISCHARGE cusecs 2DOO0O 2501000 300000 FIG. 6.11 SHIFT OF RATING CURVE FOR RIVER INDUS AT PARTAB BRIDGE o <7\ 250,000 0 30 60 90 120 150 180 210 240 270 300 330 360 TIME days FIG. 6.12 ANNUAL DISCHARGE HYDROGRAPH OF RIVER INDUS AT PARTED BRIDGE (1979 DATA) o 108 the bare surface of the glacier, inducing a phase of recession flow. As snow cover melts, an increased area of ice becomes exposed, ablation of which increases discharge. Subsequently discharge falls after a fresh snowfall covers the glacier surface again. 6.5 Link Model Results The Link Model was setup for the River Indus main stem above Tarbela dam i.e, between the gauging stations at Kachura and Besham Qila (see Fig 6.4). This segment of the river was divided into three reaches namely Kachura to Partab bridge, Partab bridge to Shaital, and Shaital to Besham Qila. There are two gauged tributaries entering into the River Indus. River Gilgit is entering into the River Indus approximately 5 km above Partab bridge, while River Astore joining River Indus approximately 5 km downstream of Partab bridge. The Link Model was run for the available flow data period of seven years i.e from 1979 to 1985. Results generated from tiie model were compared with that of the recorded flows. For different years and for different gauging stations various plots were generated. These included (i) Observed and recorded hydrographs for overall comparison of the model performance, and detecting errors in flow rating curves. (ii) Linear plots for detecting linear relationships between observed and calculated flows. In addition to this two least squares criteria were also applied to estimate the performance of the model. The coefficient of efficiency E was calculated to measure the over all performance of the model, while the coefficient of 109 determination D was calculated to evaluate the model performance after removing linear errors. The results and explanations on these results are presented below: The annual recorded and calculated discharge hydrographs of River Indus at Partab bridge are presented in Figure.6.12a. The overall routing for this period seems to be reasonable, there is some noticeable difference at the peaks , where the calculated flows are less than the recorded flows. This difference may be due to the local effect of weather which at lower elevations produced rain and subsequently effected the upstream area of flie basin by inhibiting the solar radiation, therefore resulting in a sharp drop in flows after the event. The other noticeable feature is that flows from September to November are under estimated by the Model, this may be due to tiie ground water contribution which can be calculated by taking difference of routed and calculated flows. Figures 6.13 to 6.16 are the hydrographs of the routed and recorded flows of River Indus at Partab bridge. In these hydrographs the effect of the shifting the rating curve is prominent and is progressively increasing over the years for increased flows. This station seems to be susceptible to this effect, and Figure 6.11 shows the change of rating curve in one year for the same station. Figure 6.17 and Figure 6.18 are the plots of recorded and routed flows for the period of 1984 and 1985 respectively at the next downstream gauging station Shaital. Important feature of these graphs is the indication of ungauged lateral flows entering to the reach. Figure 6.19, is a plot of calculated and recorded hydrographs at Besham Qila and it shows that in this reach large volumes of ungauged lateral flows are entering. The form of the recorded hydrograph has a remarkable similarity with that of the routed hydrograph. This shows that the ungauged lateral flows are affected by the weather pattern no in a similar way as the routed flows, and there is a similarity in the response of the basin as well. Ungauged lateral flows were calculated by taking difference of the recorded and calculated flows and were plotted with the hydrograph from the neighborhood sub-basin Astore and is shown in Figure 3.16. This figure also shows some similarity in the timing and distribution of both hydrographs. The estimated lateral flows show greater variability than the Astore data, but the general pattern is similar. In addition to the hydrograph comparison as discussed above, linear plots of calculated and observed flows were also drawn to detect linear relationship between observed and forecasted flows. Figure 6.20 is a linear plot of recorded and calculated flows of River Indus at Partab bridge for the period of 1979. This figure shows that, in general, very large flows are imder estimated, while very small ones are closely estimated and medium ones are best estimated. Coefficient of efficiency and coefficient of determination of the routed flows for different period at Partab bridge and Besham Qila gauging station were also determined. For 1979 flows the coefficient of efficiency E and coefficient of determination D at Partab bridge were 99.19% and 99.37% respectively, and at Besham Qila gauging station E and D values were 28% and 92 % respectively. This clearly shows that between Kachura and Partab bridge small volumes of ungauged lateral flows are entering and routed and recorded flows are highly correlated, whereas, the 28% E value for the routed flows at Besham Qila indicates that large volumes of ungauged lateral flows are entering between Kachura and Besham Qila, However the 92% coefficient of determination for flows between Kachura and Besham QUa indicates that these ungauged lateral flows are highly correlated with the main stream flow. As the value of D is quite close to one the difference of recorded and I l l calculated flows is a good estimate of ungauged lateral flows entering between Kachura and Besham Qila. 6.6.. Conclusions The Literature review reveals that the channel routing is essentially a solution of Saint-Venant Equations 3.1 and 3.2 for unsteady open channel flow. Depending on the nature of the channel and the flow, various simplifications and assumptions are often made to these equations to suit the conditions under study, and to reduce the complexity of the calculations. The results given in Section 6.5 show that the method proposed by Quick and Pipes (1975) for non-linear channel routing can effectively and easily be applied, with little adjustment, to the situation prevailing in the study area. This method reasonably models the routed flows and, in addition, the method can be used to estimate the ungauged lateral flows, which can create considerable problems in some routing methods. It is shown that, if the stream gauging data describing the channel flow characteristics and representing tiie whole channel is provided, the Flow Model can be used to determine these ungauged lateral flows. In the present study these flows were evaluated by applying statistical methods and by comparing these flows with the neighboring hydrologicaUy similar Astore watershed hydrograph. Application of the Flow Model to the river system as a whole was done with the help of the Link Model, which was used to connect aU the tributary watersheds inputs, channels, sub-reaches into a main channel system. At the present time, for this study, tributary watershed inputs to the Link Model are the recorded flows at the gauging stations, but in the future this model can 112 easily be linked with the UBC Watershed Model, which can be used to estimate the tributary watershed flows. Performance of the Flow Model was evaluated by plotting observed and routed hydrographs and by using several least square criteria. The coefficient of efficiency E was used to measure tfie overall performance of the Model, while evaluating the Model performance by using the coefficient of determination D gives an evaluation which is not influenced by linear errors. Therefore potentially removable errors can be estimated by taking the difference of the two criteria, D-E. The same statistical criteria have also been used to check the validation of estimated ungauged lateral flows. The present study therefore shows that within the accuracy of the available data, flow routing can be carried out with reasonable accuracy, as is demonstrated at Partab bridge when there is little ungauged lateral flows entering between Kachura and Partab bridge. Even when large ungauged lateral flows exist, between Partab bridge and Besham Qila flie model is shown to work reasonably well, and can be used to give estimates of these ungauged lateral flows, if the downstream outflows are known. These estimated ungauged lateral flows can be used to assist in calibrating the watershed model for estimation of these tributary flows. For real-time forecasting, the various upstream inflows and lateral tributary flows would need to be estimated using a watershed model, or measured at the tributary inflow points. The Flow model has been shown to be a suitable method for routing all these inflows to the downstream point, which then, for this study, becomes the estimated reservoir inflow to Tarbela Dam. 250,000 S 0 m 200,000 150,000 Jg 100,000 50,000 LEGEND RECORDED FLOWS CALCULATED FLOWS FIG. 6.12a RECORDED AND CALCULATED FLOW HYDROGRAPHS OF RIVER INDUS AT PARTAB BRIDGE u> 3001000 o w 2901000 2001000 1901000 o 1001000 501000 LEGEND RECORDED FLOWS ROUTED FLOWS 30 eO 90 120 190 180 210 240 270 300 330 360 TIME days FIG. 6.13 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT PARTAB BRIDGE (1980) 3601000 3001000 250i000 200^000 8 3 I 1S0IOOO o w D 1001000 501000 0 LEGEND RECORDED FLOWS ROUTED FLOWS 0 30 60 90 120 150 180 210 240 270 300 330 360 TIME days FIG. 6.14 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT PARTAB BRIDGE (1983) t / i Ui o 4001000 38D100O 3001000 2S01OOO 2001000 1501000 1001000 901OOO LEGEND RECORDED FLOWS ROUTED FLOWS 0 30 eO 90 120 190 180 210 240 270 300 330 360 TIME days FIG. 6.15 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT PARTAB BRIDGE (1984) o\ 2501000 2001000 1501000 8 Si 3 g 1001000 o w o SDlOOO LEGEND RECORDED FLOWS ROUTED FLOWS I I J I 30 60 90 120 150 180 210 240 270 30O 330 380 TIME days FIG. 6.16 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT PARTAB BRIDGE (1985) 300KXX) o w 25D1000 2001CXX) i 1501000 1001000 501000 r LEGEND RECORDED FLOWS ROUTED FLOWS 30 60 90 120 150 180 210 240 270 300 330 360 TIME days 00 FIG. 6.17 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT SHATIAL (1984 DATA) 2SDIO0O 8 8 3 UJ Q O « Q 20QOOO ISOiOOO 1001 OX) 501000 RECORDED FLOWS ROUTED FLOWS I J 0 30 60 90 120 190 180 210 240 270 300 330 360 TIME days ^ FIG. 6.18 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT SHATIAL (1985 DATA) 4001000 HI o w 3001000 2001000 1001000 O FIG. 6.19 COMPARISON OF ROUTED AND RECORDED HYDROGRAPHS OF RIVER INDUS AT BESHAMQILA( 1979) 250.000 200.000 150.000 o 2 Q LU i 100.000 50.000 ..-•••(5)0 O o %•% o oo o c9 d o O" --0, o .---• o P"' o . LINE OF PERFECT MATCH ± 50.000 100,000 150.000 200.000 250.000 OBSERVED FLOWS cusecs FIG. 6.20 LINEAR PLOT OF OBSERVED AND CALCULATED RIVER INDUS FLOWS AT PARTAB BRIDGE (1979 DATA) N) 122 REFERENCES Abbott, M. B. (1966) An Introduction to the Method of Characteristics. New York: American Elsevier. Bradely,C. (1990) "Snow and Ice-melt as Determinants of Hydrology in the Upper Indus Basin, Pakistan" In Snow and Ice Hyrology Project, Uppeu Indus Basin Vol 4, 1990 Cunge, J. A. (1969) " On the subject of a flood propagation computation method(Muskingum method)" J. Hydraulic Res., 7 (2), 205-230. Harley, B. M., Perkins, F. E., and Eagleson, P. S. (1970) A Modular Distributed Model of Catchment Dynamics. Hydrodynamics Lab. Report 133, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts. Harris, G. S. (1970) "Real time routing of flood hydrographs in storm sewers". J. Hydraulic Div., ASCE, 96 (HY6), 1247-1260. Hayami, S. "On the Propagation of Flood Waves". Bulletin No.l, Disaster Prevention Research Institute, Kyoto University, Japan(December 1951). Henderson, F. M. (1966) Open Channel Flow. New York: Macmillan Co., Chap 8 and 9. pp. 285-394. Hewitt, K. (1982) " Natural dams and Outburst floods of the Karakoram Himalaya". In: Hydrological Aspects of Alpine and High Mountain Areas. Glen, J.W. ed., lAHS Pubhcation No. 138, pp. 259-269. 123 Hewitt, K. (1985) Snow and ice hydrology in remote, high mountain areas: The Himalayan sources of the River Indus. Snow and Ice Hydrology Project Working Paper. No.l, Wilfrid Laurier University, 29 p. Hewitt, K, and Young,G.J (1988) "Hydrology research in the upper Indus basin, Karakoram Himalaya, Pakistan". In Hydrology of Mountainous Areas. IAHSPubl.no. 190, 1990. Liggett, J. A., and Cunge, J. A. (1975) " Numerical methods of solutions of the unsteady flow equations". In Mahmood, K., and Yeyjevich, V. (eds). Unsteady Flow in Open Channels, Vol. I. Fort Collins, Colorado: Water Resources Publications, Chapter 4, pp. 89-182. Linsley, R. K., Kohler, M. A., and Paulhus, J. L. H. (1982) Applied Hydrology. New York: McGraw Hill Book Co., pp. 266. Linsley, R. K., Kohler, M. A., and Paulhus, J. L. H. (1949) Applied Hydrology. New York: McGraw Hill Book Co., pp. 502-530. Linsley, R. K. (1971) A Critical Review of Currently Available Hydrologic Models for Analysis of Urban Storm Water Runoff. Report, Hydrocomp International Inc., Palo Alto, California. Mayewski, P.A., and Jeschke, P.A. (1979) "Himalayan and Trans-Himalayan glacier fluctuations since AD 1812"-/. Arctic and Alpine Research Vol. 11, pp. 267-287. Mayewski, P. A., Pregent, G. P., Jeschke, P.A., and Ahmad, N. (1980). "Himalayan and Trans-Himalayan glacier fluctuations and the South Asian monsoon record. / Arctic and Alpine Research. Vol. 12, No. 2, pp. 171-182. Miller, W. A., and Cunge, J.A. (1975) "Simplified equations of unsteady flow". In Mahmood, K., and Yeyjevich, V. (eds) Unsteady Flow in Open 124 Channels, Vol. I. Fort Coolins, Colorado: Water Resources Publications, Chapter 5, pp. 183-257. Nash, J.E. and Sutchffe, J.V. (1970) " River Flow Forecasting Through Conceptual Models, I. A Discussion of Principles" , Journal of Hydrology, Amsterdam, Vol. 3, pp. 282-290. Tatum, F. E. (1940) A Simplified Method of Routing Flood Flows Through Natural Valley Storage. Unpublished memo, US Engineers Office, Rock Island, Illinois. Quick, M. C , and Pipes, A. (1975) "Nonlinear channel routing by computer" . J.Hydraulic Div., ASCE, 101 (HY6), 651-664. Saint-Venant, B. De (1871) "Theory of unsteady water flow, with applications to river floods and to propagation of tides in river channels",. Computes Rendus Acad. Sci., Paris, Ti, 148-154, 237-240. ( Translated into English by US Corps of Engineers, No. 49-g, Waterways Experiment Station, Vicksburg, Mississippi, 1949.) Wake , C.P. (1987) "Spatial and temporal variation of snow accumulation in the Central Karakoram, Northern Pakistan". MA thesis. Dept. of Geography, Wilfrid Laurier University, pp. 121-122. 125 APPENDICES • *** BASIN HOOEL •*• • This model routes channel flows through different reaches and in the 'process adds lateral I flous to the main channel flows. Addition of lateral •floMS is usually carried out at the end of the reach, but if necessary they 'can be added at the begining of the reach. This Model uses three types of files ' (1) 8asin.dat: This file contains all necessary information regarding 'nodenuitiers, nunber of reaches, where to add lateral flows and the name of 'constant file. This ASCII file can be prepared from HAKBASIN.DAT. • (2) Constant.dat: This ASCII file contains reach characteristics data such •as reach length, area discharge and velocity discharge relation ships at the •nodes, storage constant modification factor,flow data time interval etc,. • (3)Station.flw : These files are flow files in binary format, and can be •prepared fron HAKEINFL.EXE. • zsi i«axKaaaaBsaniBsSBMxa«xKBaisaassasKaiBB3iBSBS3WaaBSBBaBBBSsamxnxBSSss«ssns««s«B •SOYNAMIC COMMON SHARED iS TYPE FLOWTYPE DTE AS SINGLE a AS SINGLE END TYPE DECLARE FUNCTION MESSACE1S <) {• DECLARE FUNCTION MESSAGE2S () *B DECLARE FUNCTION JulZOateS (JDAYlf) 'O DECURE FUNCTION JutianI (Yl, HI, Dl) S P a DECLARE SUB LAKEROUTING <iS, Nl, RNI, DTI, CHNLFLOW AS FLOWTYPE, LATFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE) H DECLARE SUB REACH0ATA2 (SSEGSI) ).« DECLARE SUB REACHDATA1 (RLI, RNI, ASEGSI, VSEGSI, DTI, RTKI, STKI) DECLARE SUB INSTRUCTIONS (GUIDES, CNSTFS, RN, N, iS, FJDI, RnameS, CHNLFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE, LATFLOW AS FLOWTYPE) DECLARE SUB PRINTGRAPH (xX, YX, N, QOTIO) DECLARE SUB LPRINTOUT (RN, Nl, Q U O , Q O K ) , QOTI(), TTIO, Vl(), DTEIO) DECLARE SUB SCREENOUTPUT (xX, RN, N, Q U O , OOIO, QOTIO, TTIO, Vl(), DTEO) DECLARE SUB INFLOWS (N, iS, FJDI, RnameS, CHNLFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE) DECLARE SUB ADDLATFLOW (iS, LATFS, FJOl, RN. N, CHNLFLOW AS FLOWTYPE, UTFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE) DIN SHARED CHNLFLOW AS FLOWTYPE DIM SHARED UTFLOW AS FLOWTYPE DIM SHARED IFLOW AS FLOWTYPE INPUT ••FDT^ ', FOT INPUT ••LOT^ ,^ LDT CLS OPEN '•indbasin.dat" FOR INPUT AS #30 LINE INPUT « 0 , DUMMYS 'File description INPUT #30, RnameS •River name INPUT #30, NH -Node nusbers TRCH • HN - 1 •Total reaches DECLARE SUB TIMECONVERT (YYSI, MMSI, DDSI, YYEI, MMEI, DDEI, FDTI, LDTI) DECLARE SUB CHANNELROUTING (iS, RnameS, Ni, RNI, RLI, ASEGSi, VSEGSI, DTI, RTKI, STKI, CHNLFLOW AS ANY, UTFLOW AS ANY, IFLOW AS ANY) ON Fi le : A:\M0OEL\BASNH0O1.BAS Page 1 of 16 INPUT #30, CNSTFS 'Constant file name INPUT #30, j$ 'Inflow file name CALL TIMECONVERT(YyS, HHS, DOS, YYE, HHE, ODE, FDT, LOT) FJD « Julian(YYS + 1900, HHS, DOS) LJD • Julian(YYE + 1900, W E , DOE) N » LJD - FJD • 1 'No. of flows to be routed IF N < 0 THEN PRINT MESSAGEIS END ELSE CALL INFL0WS(N, i$, FJDI, RnameS, CHNLFLOU, IFLOU) OPEN CNSTFS FOR INPUT AS #40 'LINE INPUT #40, OUMHYS INPUT #30, GUIDES SZ a (N • NN • 2) • TRCH RN • 0 DIM SHARED QO(SZ), QI(SZ • 3), OOT(SZ), OOS(SZ), OTE(SZ), 00(20 • NN), OP<20 * NN), QA(20 * NN), OAOO(20 * NN), DVDQ{20 • NN), QS(20 * NN), DSOQ(20 • NN), TT(NN • N), V(NN • N) CALL INSTRUCTIONS(GUIDES, CNSTFS, RN, N, iS, FJDI, RnameS, CHNLFLOW, IFLOW, UTFLOW) END IF END REH SSTATIC SUB AODLATFLOW (iS, UTFS, FJDI, RN, N, CHNLFLOU AS FLOWTYPE, UTFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE) <aasaaasaaaaaaaaaBBBsaBaaaBaaaasBBaBBaSBaBaaaaBaaaBaBaasasaBa3BBaaB8aaB8BBBaBaaxaaaaBBBaaaBa«aaaBa8Ba • LATERAL FLOWS WILL BE ADDED TO THE HAIN CHANEL FLOW IN THIS PROCEDURE I 'ssBESBSBasaaaaaaaasaBBnaaBasaBKeaBnasaBaBvasasaaaBaasaaaaasasBaaBBaBaaaaaaaBaBaBBBsaBaaaBaBaasaEaaasa DOS » "." + LTRIMS(STRS(RN)) FLES a "ADDLAT" • 00$ St » 1 OPEN iS FOR RANDOM AS #2 LEN > LEN(CHNLFLOW) OPEN LATFS FOR RANDOM AS #3 LEN - LEN(UTFLOW) OPEN FLES FOR RANDOM AS #4 LEN » LEN(IFLOW) GET #2, , CHNLFLOW GET #3, , LATFLOW dtel » CHNLFLOW.DTE dte2 • LATFLOW.OTE diff • dtel - dte2 REDIM QL(dfff + N) FOR It » (dfff + 1) TO (diff + N) GET #3, k, LATFLOW QL()t) ' LATFLOW.Q QI(k) > OL(k} + CHNLFLOW.Q IFLOW.OTE « dtel IFLOW.O a ai(k) PUT #4, St, IFLOW •-* GET #2, St, CHNLFLOW JT File: A:\M00EL\BASNM001.BAS Page 2 of 16 dtel « dtel + 1 St = St + 1 NEXT k CLOSE 2 CLOSE 3 CLOSE 4 i$ « FLEt END SUB SUB CHANNELROUTING <i$, RnameS, N, RN, RL(, ASEGS, VSEGS, DTI, RTKI, STKI, CHNLFLOW AS FLOWTYPE, LATFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE) i3Escsx3atts«s3ssaaxazatsa8:asESxttBSBSEasBKaaB83asa«3i[3niaasax3it>snBS3iisscsa=a=saxKSB«BssassxxaaBssssEas3ssssaBSKsanas3S3S3ssxECs==s=s=Ks»s 'THIS PROCEDURE WILL ROUTE FLOWS THROUGH A SIMPLE LINEAR WEDGE TYPE STORAGE CLS OPEN it FOR RANDOM AS #6 LEN ' LEN(CHNLFLOW) GET *6, , CHNLFLOW CRDTE# a CHNLFLOW.DTE SEEK #6, 1 FOR J • 1 TO N GET #6, , CHNLFLOW OUJ + 1) a CHNLFLOW.Q NEXT J CLOSE 6 'QKN + 1) a OI(N) DD$ a "." + LTRIMt(STR$<RN • 1)) FLEU a "CHNELOOT" + DDt OPEN FLEU FOR OUTPUT AS #10 'Openlns a file for conpUta output At a '• DATE INFLOWS 0/FLOW STRG Q TRNSUTED TT/DT VEL" PRINT #10, At kkk a 0 STORAGE: 00(1) » QI(2) FOR J a 1 TO N IF OUJ + 1) <» QA(1) THEN C a 1 IF OKJ • 1) > QA(1) THEN C • 2 SELECT CASE C CASE 1 DAOQ a DAOQ(I) GOTO CALCKT CASE 2 FOR f a 1 TO ASEGS - 1 IF OKJ • 1) > QA(i) AND OKJ + 1) < QA(i + 1) THEN DADO « DAOQ(i • 1) EXIT FOR ELSE END IF ^ NEXT i t o 0 0 Fi le: A:\M0OEL\8ASNH0O1.BAS page 3 of 16 END SELECT CALCKT: KT « RL • DADO • STK / DT QO(J + 1) » 00(J) + (OI(J + 1) - QO(J)) / (1 + KT) ' PRINT 00(J + 1), NEXT J QO(N) = QO(N - 1) TRANSLATION: Q0(1) » QI(2) 00(2) » 00(1) FOR J » 1 TO N IF 00(J • 1) < QP(1) THEN A « 1 IF QO(J + 1) >» QP(1) THEN A « 2 SELECT CASE A CASE 1 IF QO(J + 1) < OPd) THEN DVDQ « DV00(1) 00 « 00(1) GOTO CALV CASE 2 FOR i . 1 TO VSEGS - 1 IF QO(J + 1) >• OP(i) AND QO(J + 1) < 0P(i + 1) THEN DVDQ » DVDQ(i + 1) 00 « 00(i + 1) EXIT FOR ELSE END IF NEXT i END SELECT CALV: V(J) - (L0G(0O(J + D ) - LOC(QQ)) • DVDQ TT(J) » HTK • RL / (V(J) * DT) IF TT(J) < .5 THEN B • 1 IF TT(J) > .5 AND TT(J) < 1 THEN B » 2 IF TT(J) > 1 AND TT(J) < 2 THEN B » 3 IF TT(J) > 2 THEN B • 4 IF TT(J) « 1 THEN B • 5 SELECT CASE B to CASE 1 ^ File: A:\MO0EL\BASNM0O1.BAS Page 4 of 16 q1 » QI(J) q2 " (QKJ) + QI(J + D ) / 2 03 « QI(J + 1) 04 » <OI(J • 1) + OKJ + 2)) / 2 05 « OKJ + 2) IF OKJ + 1) < QA(1) THEM 0 « 10 IF OI(J + 1) > 0A(1) THEN D « 20 SELECT CASE D CASE 10 DADO ' DADOd) GOTO CALCULTKT CASE 20 FOR i • 1 TO ASEGS • 1 IF QKJ + 1) > OA(i) AND QI(J + 1) < QA(f * 1> THEN DADO > DADOd + 1] EXIT FOR ELSE END IF NEXT i END SELECT CALCULTKT: KT « RL • DADO * STK • 2 / DT 001 • q1 002 • 001 + <q2 - 001) / (1 + KT) 003 • 002 • (03 - 002) / (1 + KT) 004 » 003 + (04 - 003) / (1 • KT) 005 • 004 + (OS - 004) / (1 + KT) V(J) • (L0G(Q04) - LOG(QO)) * DVDQ TT(J) • RTK * RL • 2 / (V(J) * DT) 0OT(J + 2) • 004 + (1 - TT(J)) • (005 - 004) CASE 2 0OT(J + 2) - 00(J + 1) • (1 - TT(J)) * (00(J + 2) - Q0{J + 1)) CASE 3 0OT(J • 2) - 00(J + 1) + (1 - TT(J)) • (QO(J + 1) • 00(J)) CASE 4 0OT(J • 2) « Q0(J) + (2 - TT(J)) * (QO(J) - 00(J - 1)) CASE 5 Q0T(J + 2) • 00(J + 1) END SELECT DTE(J) » CRDTE# CRDTE# » CRDTE# + 1 OTE# • DTE(J) SDATE » VAL(Jul20ate$(0TE#)) kkk » kkk + 1 PRINT #10, SDATE; < ^ O File: A:\M0OEL\BASNHQD1.BAS Page s of 16 PRINT #10, » '•; PRINT #10, USING •• timmt»*.im »; OKJ * 1); PRINT #10, USING •• mmmm.int »; eocj + D ; If kkk « 1 THEN PRINT #10, " " ELSE PRINT #10, USING " mrnmt.m • ; QOTCJ + Z ) ; PRINT #10, USING " ##.### "; TT<J); PRINT #10, USING » ###.## "; V(J) END IF OTE# « DTE(J) PRNTDTE > VAL(JulZOate$<OTE#)> VIEW PRINT 20 TO 23 PRINT " « INPROGRESS >>» PRINT " REACH No.("; RN; ") ••; PRNTDTE 'PRINT QOT(J + 2): SLEEP NEXT J CLS 2 VIEW PRINT 1 TO 25 CLOSE 10 St • 1 DOS • "." + LTRIM»(STR$(RN)) FLEt • "FLOUOUT" • 0D$ FLQOJ • "FLOWSTRG" + DDJ FLQIS » "FLOWIN" * DD» OPEN FLQIS FOR RANDOM AS #61 LEN - LEN(UTFLOU) OPEN FLQO$ FOR RANDOM AS #62 LEN • LEN(CHNLFL0U> OPEN FLES FOR RANDOM AS #8 LEN = LEN(IFLOU) FOR k » 1 TO N IFLOW.OTE » DTE(k) + 1 IFLOU.Q ' OOT(k * 2) PUT #8, St, I FLOW LATFLOW.DTE • DTE(k) CHNLFLOW.DTE • DTE(k) LATFLOU.Q • Q K k + 1) CHNLFLOU.Q » QO(k + 1) PUT #61, St, UTFLOU PUT #62, St, CHNLFLOW St » It • 1 NEXT k CLOSE 8 CLOSE 61 CLOSE 62 i$ • FLES END SUB U> File: A:\MQ0EL\BASNM001.BAS p^a* i o SUB INFLOWS <N, il, FJOl, Rnamet, CHNLFLOW AS FLOWTYPE, IFLOU AS FLOWTYPE) • PREPARES CHANNEL FLOW WORKING FILE IssssaaaaaaaaaaasaaassaaxBasaascBaasaaaiEsasaxaB] DOS > "." + LTRIMt(STRS(1)) FLES « RnameS + D0$ OPEN iS FOR RANDOM AS «2 LEN ' LEN(CHNLFLOW) OPEN FLES FOR RANDOH AS #3 LEN ' LEN(IFLOW) GET #2, , CHNLFLOW DTE « CHNLFLOW.OTE :CBBB3Baxa3sasasaaBSBasBasBaaBsaxs DTE + 1 / LEN(CHNLFLOW) RECNOS THEN k » FJO P • FJD RECNOS • L0F(2) IF k < 1 OR N > LOCATE 13, 5 PRINT MESSAGE2S END ELSE SEEK #2, 1 FOR J » 1 TO N GET #2, k, CHNLFLOU IFLOW.Q « CHNLFLOW.a IFLOW.DTE - P PUT #3, J, IFLOW P • P + 1 k » k + 1 NEXT J iS » FLES CLOSE 2 CLOSE 3 END IF END SUB SUB INSTRUCTIONS (GUIDES, CNSTFS, RN, N, JS, FJOl, RnaoeS, CHNLFLOW AS FLOUTYPE, IFLOU AS FLOUTYPE, LATFLOU AS FLOWTYPE) • This procedure will select the procedure to preform the assigned tasks DO IF GUIDES IF GUIDES IF GUIDES IF GUIDES IF GUIDES IF GUIDES "CHANNEL" THEN X » 1 "LAKE" THEN X " 2 "ADDLATL" THEN X « 3 "NEUREACH" THEN X « 4 "FILEOUT" THEN X » 5 "JOBENO" THEN X » 6 'To perform channel routing 'To perform take routing 'To add lateral flows 'To assign new reach nuB^r 'To prepare output file 'To exit from the program SELECT CASE X CASE 1 CALL REACHDATAKRL, RN, ASEGS, VSEGS, DT, RTKI, STKI) File: A:\M00EL\BASNM001.BAS Page 7 of 16 CALL CHAHMELROUTING(iJ, RnameJ, N, RN, RLI, ASEGS, VSECS, DTI, RTKI, STK!, CHMLFLOU, LATFLOW, IFLOU) INPUT #30, GUIDES IF LTRIMSCGUIDEt) • "JOBENO" THEN EXIT DO END IF CASE 2 CALL REACHDATA2(SSEGS) CALL LAKEROUTlNG(i$, N, RN, DTI, CHNLFLOW, LATFLOW, IFLOW) INPUT #30, GUIDES IF LTRIMS(GUIDE$) ' "J08END" THEN EXIT DO END IF CASE 3 INPUT #30, LATFS CALL ADOLATFLOWdt, LATFJ, FJOl, RN, N, CHNLFLOW, LATFLOW, IFLOW) INPUT #30, GUIDES IF LTRIHJ(GUIDEI) • "JOBEND" THEN EXIT DO END IF CASE 4 RN ' RN + 1 INPUT #30, GUIDES IF LTRIMJ<GUIDEJ) • "JOBEND" THEN EXIT DO END IF CASE 5 NN > RN + 1 DOS « "." + LTRIMKSTRJCNN)) FOUTS ° RnameJ * DDS OPEN (J FOR RANDOM AS #S0 LEN • LEN(IFLOU) OPEN FOUTJ FOR ftAHDOH AS #60 LEN ' LEN( LATFLOW) FOR i « 1 TO EOF(50) GET #50, J, IFLOW LATFLOU.DTE - IFLOW.DTE LATFLOW.Q « IFLOW.O PUT #60, I, LATFLOW NEXT i CLOSE 50 CLOSE 60 INPUT #30, GUIDES IF LTRIMS(GUIOEI) « "JOBEND" THEN EXIT DO END IF CASE 6 IF LTRIM»(GUIOEJ) • "JOBEND" THEN EXIT DO END IF END SELECT LOOP File: A:\MO0EL\BASNH0O1.BAS Page 8 of 14 END sua FUNCTrON JutZDateS (JDAY#) ' This function returns YYMHDD date f o m t from Julian date •xa9BXB3aa3Ba333xssB3s833as3aaaBVBa3axtt«aaaasxa3a3a3S3c«38aass8»asaB3asaaxa3a Start • FIX(JDAY# / 365.25) 'compute a start year to save time FOR Y » start TO start + 100 CJDay 3 (365 • (Y - 1)) * FIX((Y - 1) / 4) - FIX((Y - 1) / 100) • FIX((Y - 1) / 400) + 1 IF CJDay > JDAY# THEN YEAR » Y - 1 EXIT FOR 'the correct year has been found Hhen cjday>jday END IF NEXT Y Y B YEAR 'the correct year •codfjute leap year correction IF (Y / 4 a FIX(Y / 4) AND NOT (Y / 100 » FIX(Y / 100))) OR Y / 400 a FIX(Y / 400) THEN C • 1 ELSE C a 0 END IF •coirpute day of year YOAY • JDAY# - ((365 • (Y - 1)) + FIX((Y - 1) / 4) - FIX((Y - 1) / 100) • FIX((Y - 1) / 400)) 'coopute month ninber FOR H B 1 TO 13 IF H < 3 THEN CYDAY » CINT((M - 1) • 30.58333) + 1 ELSE CYDAY » CINT((M - 1) * 30.58333) - 2 + C • 1 END IF IF CYDAY > YOAY THEM MONTH B M - 1 EXIT FOR 'the correct month has been found when cyday>yday END IF NEXT M H • MONTH 'the correct month •compute day of the month IF M < 3 THEN D • YOAY - CIMT((M - 1) • 30.5833) ELSE D « YOAY - (CINr((H - 1) * 30.58333) - 2 • C) END IF ystrS 3 LTRIM*(STR$(Y MOO 100)) IF M < 10 THEN monstr* ' "0" * LTRIM*(STR$(M)) ELSE monstr$ • LTRIH$(STR$(M)) END IF IF D < 10 THEN daystrS » "0" + LTRIH$(STR$(D)) ELSE ^ daystrS ' LTRIMS(STRS(D)) U> ••> File: A:\M00EL\BASNM0D1.BAS p,g END IF Jut20ate$ " ystr$ • monstrS + daystrS END FUNCTION FUNCTION Julian (Y, M, D) • CONVERTS YYWfflD DATE FORMAT TO JULIAN DATE •coinpute the leap year correction IF (Y / 4 » FIX(Y / 4) AMD MOT (Y / 100 • FIX(Y / 100))) OR Y / 400 • FIX(Y / 400) THEN C « 1 ELSE C • 0 END IF 'coofXite the day of the year IF M < 3 THEN YDAY • CIMT<(M - 1) • 30.58333) + D ELSE YDAY • CINT((M - 1) * 30.58333) - 2 + D + C END IF •conpute the Julian day Julian • (Y - 1) • 365 + FIX((Y - 1) / 4) - FIX((Y - 1) / 100) + FIX((Y - 1) / 400) + YDAY END FUNCTION SUB LAKEROUTING (i$, N, RM, DTI, CHNLFLOW AS FLOWTYPE, UTFLOW AS FLOWTYPE, IFLOW AS FLOWTYPE) • THIS PROCEDURE WILL ROUTE THE CHANNEL FLOW THROUGH A LAKE CLS REDIH QI(N • RN), QOS(N + RH), DTE(M + RN) OPEN is FOR RANDOM AS #20 LEN « LEN(CHNLFLOW) GET #20, , CHNLFLOW CRDTE# • CHHLFLOW.DTE SEEK #20, 1 FOR J a 1 TO N GET #20, , CHNLFLOW QI(J) « CHNLFLOW.Q NEXT J CLOSE 20 QH2) = QUI ) QOS(I) « 01(1) FOR J » 1 TO N IF QKJ + 1) < QS(1) THEN C « 1 IF QKJ + 1) > QS<1) THEN C - 2 SELECT CASE C CASE 1 •-» OSOQ « DSOO(I) ^ File: A:\HO0EL\BASNN0O1.BAS Page 10 of 16 GOTO CALCKT2 CASE 2 FOR i a 1 TO M - 1 IF QI<J + 1) > OS(i) AND QI(J + 1) < QS(i + 1) THEN DSOO « DSOQ(i + 1) EXIT FOR ELSE END IF NEXT i END SELECT CALCICT2: KS ' DSOO / OT QOS(J + 1) « OOS(J) + (OI(J • 1) - OOS(J>) / (1 + KS) NEXT J CRDTE# « CROTE# + 1 DOS » "." • LTRIH$(STR$(RN)) FLEt • RnameS * DOS OPEN FLES FOR RANDOM AS «8 LEN > LEN(IFLOU) FOR k « 1 TO N IFLOW.DTE < CRDTE# IFLOW.Q B QOS(k * 1) PUT #8, k, IFLOW CRDTE# - CRDTE# * 1 NEXT k CLOSE 8 i$ « FLE* END SUB SUB LPRINTOUT (RN, N, Q U O , QOI(), OOTK), TTK), V K ) , DTEO) laaBsasaBKESBaBaBaBBsaBaBSBaasssaaatfBSBaKBaBBaaxaaaasanaasaxaKStBamBsaaaaaBaiasaa ' THIS WILL PRINT OUTPUT ON PAPER laBaBaSBaaBaaaaBBBSaaaaBaaaBBBaaBaBBBBaBaBBaBaBaBaBaaBaBaaaBBBsaBxaBaaBBBaaaaBB LPRINT " RIVER INDUS FLOW ROUTING FOR REACH No. ("; RN; ")» LPRINT LPRINT " DATE INFLOWS 0/FLOW STRG 0 TRNSLATED TT/DT VEL" LPRINT kkk = 0 CNT » 0 FOR J B 1 CNT a TO CN1 IF CNT » CNT » LPRINT " LPRINT LPRINT " LPRINT ELSE 0 N • + 1 57 THEN DATE RIVER INDUS FLOW ROUTING FOR REACH No. (••; RN; ")•• INFLOWS 0/FLOW STRG Q TRNSLATED TT/DT VEL" DTE# = DTE(J) SOATE « VAL(Jul2D«te$(DTE#)) kkk » kkk + 1 LPRINT " "; LPRINT SDATE; LPRINT " "; LPRINT USING •• mmim.m "; QI(J); f^ o\ File: A:\MODEL\BASNM001.8AS Page 11 of 16 LPRINT USING " immm.itit ••; QO(J); IF kkk • 1 THEN LPRINT •• " ELSE LPRINT USING " mmmi.m "; QOT(J + 1 ) ; LPRINT USING " ##.##« "; TT(J); LPRINT USING " ###.## "; V(J) END IF END IF NEXT J END SUB FUNCTION HESSAGEU 'BssaaBasBaBBaBaaaaBaaaaaBaaBaasaaxBaaaaaBBBsasaBBBaaaaBaaaaaaaBxaBaaBaxaaaaaaB ' This wUl print the following message •KsaBaaaBanaBBaaBasasBBaaaassaaasaaaaaaaaBBasBaaasBBaaaaBaasBBaaxaBaEBaassBaBax CLS LOCATE 13, 5 MESSAGE1S « » STARTING DATE SHOULD COHE FIRST" END FUNCTION FUNCTICM MESSAGE2$ insaaaaaaaaxaaaBaBBaaaaBBaaaaaaBaaaBSBaBBasaBaBaaBaaBxasBBBBaasaasaaaasaxaBaaBa ' This will print the following message •axaEaxEBBsaaaaaaasaaBaBBBaaBaansBBaasaaaaBaaaaaiucasaBaaaaaaaBBaBxaBBxasxBaaBBB CLS LOCATE 13, 5 MESSAGE2$ » "ROUTING DATES ARE NOT MATCHING INFLOW FILE DATES" END FUNCTION SUB PRINTGRAPH (xX, YX, N, QOTK)) 'This procedure Hflt produce a hydrograph froB values calculated fn the last reach <aBSBBassBXBaaBaaxxXBaBaaaasaEaB=SBBaaBaaBaBaaBaaaB3aaaBaaaaaaBBaaaa«asa3B8B8aa CLS LOCATE 13, 3 PRINT "Do you want to see graphic output fro« Reach No.("; xX; ") Y/N" PRINT INPUT » "; GRPHS GRPHS » UCASE$<GRPH$) IF GRPH$ » "T" THEN DIM YAXISCIS), R<15), XAXIS(20), C(20) DIM x(2 • N + 1), Y(M + 1 ) ,_* DIM THAX(M + RN), TMIN(H • 1), YMAX<M + RN), YMIH<N + 1 ) t>> -J File: A:\HOOEL\BASNMOD1.SAS Page 12 of 16 TMAX » H: TMIN - 1: YMAX « -10: YMIM « 9999999 FOR i => 1 TO N x(i) « i NEXT f FOR i • 1 TO N x(i) » INT<x(f)) Y(i) « INT<QOT<i)) NEXT i FOR i " 1 TO N IF Y(i) >• YNAX THEN YMAX » Y(i) IF Y(i) <» YMIN THEN YMIN » Y(i) NEXT i FOR J « 1 TO 11 YAXIS(J) • YMIN + INT((YMAX - YMIN) / 11 * J) NEXT J FOR J « 0 TO 12 XAXIS(J) « (INT<((TMAX - TMIN) + 1) / 12) • J) NEXT J FOR i • 1 TO N x(i) . 80 + INT(510 • (x(f) - TMIN) / (TMAX - THIN)) Y(i) « 408 - INT(352 * (Y(f) - YMIN) / (YMAX - YMIN)) NEXT f SCREEN 12 LINE (80, 56)-(590, 56), 5 LINE (80, 408)-(590, 408), 5 LOCATE 2, 12: PRINT ,- CAPTIONS •TOP LINE 'SPECIFIES X-AXIS FOB IX « 37 TO 590 STEP 42 LINE (IX, 408)-(IX, 403) NEXT IX LOCATE 2, 10: PRINT " LOCATE 30, 10: PRINT " k • 10 FOR i » 0 TO 12 C(i) • k • f * 5 LOCATE 27, C(i) PRINT XAXIS(i) NEXT f LINE (590, 56)-(590, 408), 5 LINE (80, 56)-(80, 408), 5 FOR lY » 56 TO 408 STEP 32 LINE (80, IY)-(85, lY) NEXT lY LOCATE 7, 1: PRINT "D" LOCATE 8, 1: PRINT "I" ROUTED FLOWS OF RIVER INDUS FROM REACH No. ("; xX; T I M E ( days )"; •LEFT MARGIN 'SPECIFIES Y-AXIS 'PRINTING Y-AXIS LABEL Fits: A:\M0DEL\BASNM001.BAS Page 13 of 16 LOCATE 9, 1: LOCATE 10, 1: LOCATE 11, 1: LOCATE 12, 1: LOCATE 13, 1: LOCATE 14, 1: LOCATE 15, 1: PRINT "S" PRINT "C" PRINT "H" PRINT "A" PRINT "R» PRINT «a« PRINT "E» LOCATE 17, 1: PRINT "c» LOCATE 18, 1: PRINT "f" LOCATE 19, 1: PRINT "s" LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE LOCATE 4, 6, 8, 10, 12, H, 16, 18, 20, 22, 24, 26, 3: 3: 3: 3: 3: 3: 3: 3: 3: 3: 3i 3; PRINT YMAX PRINT YAXIS(IO) PRINT YAXIS(9) PRINT YAXIS(8) PRINT YAXIS(7) PRINT YAXIS(6) PRINT YAXIS(5) PRINT YAXIS(4) PRINT YAXIS(3) PRINT YAXIS(2) : PRINT YAXIS(I) ; PRINT YHIN 'PRINTING Y-AXIS VALUES FOR i » 1 TO N - 1 LINE (x(i), Y(i))-(x(f + 1), Y(1 • D ) , U , , IHBBBB NEXT i • LINE (80, 470)-(1S8, 470), 14, , tHBBBB DO LOOP UNTIL INKEYS <> •"• • PRINTING GRAPH ELSE PRINT END IF SCREEN 0 END sua SUB REACHDATA1 (RL, RN, ASEGS, VSEGS, DT, RTK, STK) ' This procedure wUl read the reach characteristic data froa constant file LINE INPUT #40, DUMHYS LINE INPUT #40, DUMHYS INPUT #40, RL, VSEGS, ASEGS, DT, RTK, STK 'PRINT ASEGS, VSEGS: SLEEP: END REDIM aA(ASEGS), DADQ(ASEGS) LINE INPUT #40, DUHHYt REDIH 0<3(VSEGS), QP(VSEGS), DVDa(VSEGS) FOR i » 1 TO VSEGS INPUT #40, QQ(i), QP<i), DVDQ(i) NEXT i LINE INPUT #40, DUHMYS Fite: A:\M0OEL\BASNM0D1.BAS Page 14 of 16 FOR j B 1 TO ASEGS INPUT #40, QA(i), DADQ(i) NEXT i END SUB SUB REACHDATA2 (SSEGS) This procedure wjU read lake constant data LINE INPUT #40, DUMMYI LINE INPUT #40, DUWYS INPUT #40, SSEGS FOR i a 1 TO SSEGS INPUT #40, aS(i>, 0S0Q(O NEXT f END SUB SUB SCREENOUTPUT (xX, RN, N, OIIO, OOIO, QOTK), TTIO, VIO, DTEO) •KsaBaxaasaBaaaaaaaasaaaaBBBBaSBaaaBBaaaaBBaaaaaBaaaaasaaaaasaaaaBaBaaBai ;saaaaaaaa It will produce tabular output on screen •sssaaaaaBEEaasBaaaaaaaBaaaBaasaaEaaaBBaaBaBsEaaBBaa BEaBasaBSEaEaasEBaaaBBas CLS LOCATE 13, 3 PRINT "Do you want to see out put on screen froB Reach No. C ; xX; ") PRINT INPUT " SCREENS > UCASE$(SCR£ENS) IF SCREENS B Myii JHEN A « 1 IF SCREENS <> "Y" THEN A • 2 SELECT CASE A CASE 1 CLS CC a 0 kklc » 0 Y/N» "; SCREENS PRINT " DATE INFLOWS PRINT FOR J « 1 TO N IF CC • 20 THEN B « 10 IF CC <> 20 THEN B « 20 SELECT CASE B CASE 10 PRINT PRINT " *•" 0/FLOW STRG 0 TRNSLATED TT/DT VEL" PRESS AMY KEY TO CONTINUE •*•*" DO LOOP UNTIL INKEYS <> »'< CLS PRINT " DATE PRINT INFLOWS O/FLOW STRG Q TRNSLATED TT/DT VEL" o File: A:\M00EL\BASNH001.BAS Page 15 of 16 CC = 0 CASE 20 OTE# • DTE(J) SOATE = VAL(Jul20«teJ(0TE#)) kkk « kkk + 1 PRINT SOATE; PRINT " "; PRINT USING " mmmm.lHt »; OUJ); PRINT USING » mmnm.m »; oo( j ) ; IF kkk » 1 THEN PRINT " " ELSE PRINT USING •• mtmm.t* »; OOTCJ + i> ; PRINT USING • ##.### "; TT(J); PRINT USING " ###.## <•; V(J) END IF CC » CC + 1 END SELECT NEXT J PRINT PRINT INPUT " PRESS P FOR PRINTER OR PRESS ENTER TO CONTINUE"; PJ PS » UCASE»(P$) IF P$ « "P" THEN CALL LPRINTOUT<RN, N, QIK), QOK), OOTIO, TTI{), Vl(), DTEO) ELSE END IF CASE 2 PRINT END SELECT END SUB SUB TIMECONVERT (YYS, WS, DOS, YYE, MHE, DDE, FOT, LOT) ' This procedure cofwerts YYMMDD date format to year, month and days >3BS«ssx3B3B3BBaBaaBSB3a3aaB3SB3Ba3B3aaaaaBaBSKaaaaBa3BsaBsaaBBSBa3BBSaBaB8BaaK YYS « INT(<FDT) / 10000) YYE " INT((LDT) / 10000) HHS = (INT{(FDT) / 100)) • MME « <INI((LDT) / 100)) -DOS « FDT - (YYS * 10000) DOE » LOT - (YYE • 10000) END SUB (YYS * 100) (YYE * 100) - (MHS * 100) - (KME * 100) File: A:\mDEL\BASNH001.BAS Page 16 of M 142 Appendix 2 Description of Basin.dat File Following is the description of the Basin.dat file prepared for the main stem Upper Indus River and given the name as Indus.dat. Figure 4.1. is an output of this file and was prepared with the help of schematic diagram of the area under study which is shown in Figure 4.2: First five lines of the Indus.dat file are common to all basin files . First line of the file is die description of the file and is only for the reader use. Second line describes the name of the channel, in this case River Indus, this name is used for all the main output files generated from the system operation. Output files are identified by the extension, which is given to them by the node number at die downstream of the channel. Third line describes the number of nodes, this is used by the program to calculate number of reaches in the system. Fourth line describes the name of die constant fQe which includes aU the necessary constant data used in channel or lake routing. This data is arranged separately for each reach. While preparing this file care is exercised tiiat the sequence of constant data for each reach should have a direct correspondence with reach sequence i.e, first reach data at the top and the last reach data at bottom. Fifth line describes the name of the file which contains the historic inflows to the chaimel at node 1. This file is opened and a working file is prepared by extracting the inflows and the corresponding dates for these inflows, for a period prescribed by the user in the beginning. All rest of the hnes in Basin.dat file are basically instructions to the main program. In the case of Indus.dat the next instruction is Newreach, the program 143 will update the value of variable reach number (RN) to 1, (by default RN=0). The next instruction is "channel" , when this instruction is read by the program it will call a sub-program "channel". This sub-program will open the constant file , in this case INDSCNST.DAT file, and will read all the values required for channel routing for this reach. It will first perform routing of the inflows through the linear storage, and then translate these flows through the channel as described by the Flow Model in Chapter 3. In order to see the effect of the routing, an ASCII file COMPOPUT.l is generated at this stage, this file shows inflows, storage outflows, translated flows, flow travel time and flow velocity at different days. Figure A. 1 is a part of an output fi-om the same file. In addition to this three other binary files are generated, mainly carrying the flow information and are required by the graphics and statistics modules of the Link Model. Next instruction to the program is "addlatl" and the name of the file carrying these lateral flows is given on the next line of the Indus.dat file, at this point a subprogram "addlatl" opens the lateral flow file and the translated flow file prepared above. This sub-program first searches for the matching dates corresponding to the translated flows in the lateral flow file and then add these flows to the corresponding translated flows. These flows are saved in a working file. Next instruction to the program is file out, this calls a sub-program which copies the latest working file and name it as channel name.node number, in our case INDUS.2. In this manner flows are routed through different reaches in the river network until all the instructions on the basin, dat file are followed. In this process different files are generated, these files are fiulher used for different type of outputs which can be obtained by using Graphics or Statistics options, which are described in Chapter 4. 144 DATE 790101 790102 790103 790104 790105 790106 790107 790108 790109 790110 790111 790112 790113 790114 790115 790116 790117 790118 790119 790120 790121 790122 790123 790124 790125 790126 790127 790128 790129 790130 INFLOWS 8476.80 8406.16 8123.60 8017.64 8017.64 8017.64 8123.60 8123.60 8123.60 8123.60 8017.64 8017.64 8017.64 8017.64 7346.56 7947.00 7947.00 7876.36 7841.04 7805.72 7735.08 7805.72 7593.80 7523.16 7523.16 7346.56 7240.60 7099.32 6958.04 7028.68 O/FLOW STRG 8476.80 8426.73 8211.86 8074.19 8034.11 8022.43 8094.14 8115.02 8121.10 8122.87 8048.28 8026.56 8020.24 8018.40 7542.17 7829.13 7912.68 7886.94 7854.40 7819.89 7759.77 7792.34 7651.61 7560.56 7534.05 7401.15 7287.35 7154.07 7015.12 7015.12 Q TRNSLATED 8390.63 8189.74 8067.95 8032.30 8033.47 8097.41 8115.98 8121.38 8111.13 8044.92 8025.59 8019.95 7945.19 7581.04 7841.38 7908.82 7882.09 7849.31 7811.10 7764.47 7771.90 7638.88 7556.95 7516.09 7386.40 7270.71 7137.53 7015.12 7016.52 TT/DT 0.832 0.839 0.844 0.846 0.846 0.844 0.843 0.843 0.843 0.845 0.846 0.846 0.846 0.865 0.853 0.850 0.851 0.852 0.854 0.856 0.855 0.860 0.864 0.865 0.870 0.875 0.881 0.887 0.887 VEL 3.91 3.88 3.86 3.85 3.85 3.86 3.86 3.87 3.87 3.85 3.85 3.85 3.85 3.77 3.82 3.83 3.83 3.82 3.82 3.80 3.81 3.79 3.77 3.77 3.74 3.72 3.70 3.67 3.67 PIG A.l A PART OF A STANDARD TABULAR OUTPUT FROM FLOW MODEL
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Flow routing model for upper Indus river (Pakistan)
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Flow routing model for upper Indus river (Pakistan) Hashmi, Danial 1994
pdf
Page Metadata
Item Metadata
Title | Flow routing model for upper Indus river (Pakistan) |
Creator |
Hashmi, Danial |
Date Issued | 1994 |
Description | For the flow forecasting of the Upper Indus Basin (UIB) Pakistan, it is observed that large quantities of ungauged lateral flows enter many sections of the Indus River system during the monsoon and the snow melt season. The river is steep, flow travel times are short, and the data are currently available only on a mean daily flow basis. There are very few gauging stations and the travel time between them is less than a day, so that standard routing procedures such as the Muskingum method was difficult to apply because of the large amount of ungauged lateral flows entering the main channel. Calibrating the Muskingum model was therefore difficult because of the ungauged lateral flows, and also because the coefficients used in the Muskingum method are a function of routing period and vary considerably when the routing period is less or greater than the travel time. There was a need for a model which could not only route the flows through the Upper Indus reaches but also could link all the sub basins contributing to the Indus flows. The whole UIB is divided into reaches, and a digital computer Link Model is developed which links these reaches and routes the flows through these reaches using the flow routing method developed by Quick and Pipes. The Link Model adds the known lateral flows to the main channel flows and also calculates the ungauged lateral flows if the stream gauging data for the downstream station is known. This Model can easily be interfaced with the UBC Watershed Model, which forecasts the land phase of the flow and the total system of models can be used for flow forecasting of UIB. The routing method adopted is very flexible and easy to calibrate, because the principal feature of the routing method is that the routing coefficients are determined directly from the stage-discharge and area-discharge measurements at the gauging stations, so the method is pre-calibrated before commencing routing, and therefore avoids the problem of accounting for ungauged lateral inflows. The model has been tested on a 400 km segment of the main stem Upper Indus River, which lies above Tarbela dam i.e, between Besham Qila and Kachura gauging stations Fig (5.1). and was found satisfactory. Given the required data the Model can be used in any similar system. |
Extent | 4634897 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-02-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050412 |
URI | http://hdl.handle.net/2429/4943 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-0114.pdf [ 4.42MB ]
- Metadata
- JSON: 831-1.0050412.json
- JSON-LD: 831-1.0050412-ld.json
- RDF/XML (Pretty): 831-1.0050412-rdf.xml
- RDF/JSON: 831-1.0050412-rdf.json
- Turtle: 831-1.0050412-turtle.txt
- N-Triples: 831-1.0050412-rdf-ntriples.txt
- Original Record: 831-1.0050412-source.json
- Full Text
- 831-1.0050412-fulltext.txt
- Citation
- 831-1.0050412.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050412/manifest