@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Civil Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Hashmi, Danial"@en ; dcterms:issued "2009-02-23T21:16:18Z"@en, "1994"@en ; vivo:relatedDegree "Master of Applied Science - MASc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms: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."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/4943?expand=metadata"@en ; dcterms:extent "4634897 bytes"@en ; dc:format "application/pdf"@en ; skos:note "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 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 < 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 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 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$ 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 » 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 + VAL(JulZOate$ 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 \".\" + 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 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 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 > -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• 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, 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 LPRINTOUT3BS«ssx3B3B3BBaBaaBSB3a3aaB3SB3Ba3B3aaaaaBaBSKaaaaBa3BsaBsaaBBSBa3BBSaBaB8BaaK YYS « INT(