
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/UF00101079/00001
Material Information
 Title:
 A new perspective of nonlinear routing in hydrology
 Creator:
 Smith, George Francis, 1951
 Copyright Date:
 1983
 Language:
 English
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright George Francis Smith. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 821012020 ( OCLC )
ACL1882 ( LTUF ) 0030349410 ( ALEPH )

Downloads 
This item has the following downloads:

Full Text 
A NEW PERSPECTIVE OF
NONLINEAR ROUTING IN HYDROLOGY
BY
GEORGE FRANCIS SMITH
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1983
ACKNOWLEDGMENTS
This study was made possible by the National Weather Service Fulltime University Assignment Program, which funded the author's studies at the University of Florida for the 19781979 academic year.
The Director of the NWS Hydrologic Research Laboratory, a position formerly held by Dr. Eugene Peck and more recently by Dr. Michael Hudlow, has provided initial encouragement and continuing support of the work presented here.
Dr. John Schaake, Chief of the NWS Hydrologic Services Division, has guided the author's research during the past several years. His support and always carefully considered advice on both the technical and personal aspects of an effort of this magnitude are gratefully appreciated.
Dr. Wayne Huber, Professor of Environmental Engineering Sciences, University of Florida, has performed masterfully the difficult task of committee chairman for absentee doctoral research. His thoughtful review of this material and concern for both the overall technical aspects of the work and the details of presentation have substantially improved the final dissertation.
The comments and exchange of ideas provided by Dr. Konstantine Ceorgakakos, NOAANRC Research Associate, have been invaluable to the author's continuing research.
The author's colleagues in the NWS Hydrologic Research Laboratory, especially Mr. Gerald Day, have provided encouragement and stimulating discussions throughout this endeavor. Their continued expressions of concern are deserving of acknowledgment.
The careful reviews and editorial comments of Ms. Ellen Stockdale and Melanie Quinn have improved immensely the quality of this report.
The computer simulations and figures were generated on the NWS Office of Hydrology PRIME 750 computer. The author thanks Mr. Gary Kemp, of the NWS Hydrologic Services Division, for his assistance in maintaining the computer support needed for this work.
For their care and diligence with the final computer plotting of the figures, the author is especially grateful to Ms. Cynthia Brown, of the NWS Hydrologic Services Division, and to Mr. Stephen Ambrose and all the technical support staff of the NWS Hydrologic Research Laboratory.
The help of Ms. Janice Lewis, of the NWS Hydrologic Research
Laboratory, in the use of the NWS Dynamic Wave Operational Program is greatly appreciated.
For their continued confidence and encouragement, the author is grateful to his family and friends. For her efficient and efficacious assistance with the technical aspects of producing this dissertation, the author is forever indebted to his wife, Margaret. Her patience, understanding and sustaining support have been salutary to making this work possible.
TABLE OF CONTENTS
Page
ACKNOWLEDGMNTS .................. ...................... *. ii
LIST OF TABLES ...................................... *....... vii
LIST OF FIGURES ............................ ......... ... .... viii
CHAPTER
I ~A NEW PERSPECTIVE OF NONLINEAR ROUTING............. 1
Statement of the Problem.................... 1
Statement of the New Perspective.... 0....... 3
Development of the NewPerspective ...................5
SelectnonaPriarMer.................5
Calibration of the Primary Model to the Hydrologic System ..........o..o..... 6
Structure of the State Space Model .......... 6
Calibration of the State Space Model to the Primary Model ......o.......... ...... 7
The State Space Model as Surrogate for the Primary Model....... o........ o.... 8
An Application of the New Perspective .................9
Organization of the Dissertation.... .............. 11
II INSIGHTS INTO NONLINEAR ROUTING FROM
LINEAR SYSTEMS THEORY .................... ............13
Linear Systems Theory ....................... ..... ....13
Introduction......... ............................ 13
Properties of aSystem........................... 15
A General Linear Differential Equation ............19
Transformation of Moments by a Stationary Linear System ......................22
Summary................ ..................... ......25
Linearization of the SaintVenant Equations
for a General Prismatic Channel ...................26
Insights from Linear Systems Theory.. ................33
Introduction...................................... 33
The Inverse Gaussian PDF as an Analytical Solution of the Unsteady Flow Equations ........33 Moments of the Inverse Gaussian PDF .............o 34
The Inverse Gaussian PDF as an Inflow Hydrograph.......... ....... ............ 37
Additional Properties of the Inverse Gaussian PDF........................... 40
iv
A Linear State Space Approach ...................... 41
Introduction .................................... 41
Criteria for Matching Impulse
Response Functions ........................... 43
Computing the Parameters of a Gamma PDF ......... 46
A Cascade of Linear Reservoirs
State Space Model ........................... 52
Summary ................................. ....... 57
III LIMITATIONS OF THE LINEAR THEORY ...................... 58
Introduction ................... ...... ............ 58
Verification of the Linear
SaintVenant Equations .......................... 59
Verfication of the Linear State Space Model
for Small Inflow Transients ..................... 82
Comparison with the 3Parameter Gamma PDF ....... 82 Comparison with the DWOPER Model ................ 86
The Importance of Nonlinearities ................... 98
Summary ............................................ 117
IV A NONLINEAR STATE SPACE ROUTING MODEL ................. 123
Analysis of the Important Nonlinearities ........... 123
Nonlinear Characteristics of the Primary Model ..... 125
Introduction .................................... 125
Theoretical Variations with Flow Level of the Impulse Response Function for the Linearized SaintVenant Equations for Incremental Inflow Transients ............ 125
Simulated Variations with Flow Level of the Impulse Response Function for the Full SaintVenant Equations for Incremental Inflow Transients ............ 130
Properties of the Nonlinear Storage Reservoir
as a Building Block for the State Space Model... 137
Development of the Nonlinear Model ................. 141
Introduction ..................................... 141
The Final Configuration of the
Nonlinear Model .............................. 148
The Path to the Final Model ..................... 153
Summary .............................................. 158
V LIMITATIONS OF THE NONLINEAR STATE SPACE APPROACH ..... 160
Introduction ........................................ 160
Verification of the Nonlinear State Space Model.... 162
Introduction ..................................... 162
SinglePeaked Inflow Hydrographs ................ 162
A DoublePeaked Inflow Hydrograph ............... 182
Continuity .................................... 182
The Number of Reservoirs ........................ 189
CPU Time ........................................ 189
Summary ............................................ 192
VI HOW TO EMULATE ANY PRIMARY MODEL ..........*............ 193
Introduction. . ... .. so ... o .. .. ................. 193
Select a PriyMder....... e............ 193
Calibrate the PriyModl.....odel........... 193
Determine the Impulse Response Properties
of the Primary Model ........ .. .. .. ......... 195
Estimate a, n and k as Functions of Flow Level.. 197 Parameters of the Nonlinear State Space Mode l ... 197
VII RELATED PREVIOUS WORK. .......... ........ .. . .. ... 200
Related Routing Models.............................. 200
Observed Variations with Flow Level
of the System Time Lag...... .................... 203
Comparison with the MuskingumCunge Model .......... 203
VIII CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK ....... 208
Conclusions ......... .. 0 .. .. .. .. .... .. ... .. .. 208
Recommendations for Future Work........ ..... 211
APPENDICES
A A REVIEW OF SOME HYDROLOGIC ROUTING MODELS............. 216
B A FORTRAN PROGRAM TO GENERATE HYDROGRAPHS WITH
INVERSE GAUSSIAN PDF SHAPES
IN A PRISMATIC CHANNEL............ . .. ........ 224
C A FORTRAN PROGRAM FOR THE NONLINEAR STATE SPACE
ROUTING MODEL...........................227
BIBLIOGRAPHY ........................................... 233
BIOGRAPHICAL SKETCH .. ............................. .. .. 239
LIST OF TABLES
Table Title Page
3.1 Organization of Figures 3.1 through 3.18,
Comparing the Linearized and Full
SaintVenant Equations ................................ 61
3.2 CPU Times for the Numerical Solution of the
SaintVenant Equations on a PRIME 750 Computer
for Small Inflow Transients ........................... 83
3.3 Pairs of Figures with Identical Input Data
Comparing the Linearized SaintVenant Equations and the Linear State Space Model with the Full
SaintVenant Equations ................................ 90
3.4 Values of Parameters a, n and k
for the State Space Model ............................. 91
3.5 CPU Times for the Linear State Space and
DWOPER Models on a PRIME 750 Computer for
Results of Figures 3.22 through 3.27 .................. 100
3.6 Organization of Figures 3.28 through 3.43,
Comparing the Linear State Space Model
with the Full SaintVenant Equations
for Large Inflow Transients ........................... 102
5.1 Organization of Figures 5.2 through 5.17,
Comparing the Nonlinear State Space Model
with the Full SaintVenant Equations
for Large Inflow Transients ........................ 163
5.2 Pairs of Figures with Identical Input Data
Comparing the Linear and Nonlinear State Space
Models with the Full SaintVenant Equations
for Large Inflow Transients ............................ 164
5.3 Continuity Check for the Nonlinear State
Space and DWOPER Models for Results of
Figures 5.2 through 5.17 .............................. 188
5.4 The Average Number of Reservoirs Used
by the Nonlinear State Space Model
in Figures 5.2 through 5.17.......................... 190
5.5 CPU Times for the Nonlinear State Space and
DWOPER Models on a PRIME 750 Computer for
Results of Figures 5.2 through 5.17 ...................191
7.1 Values of m for Several Catchments ...................... 205
LIST OF FIGURES
Figure Title Page
1.1 An Overview of the New Perspective
of Routing in Hydrology ............................... 4
1.2 An Overview of the Development
of the New Perspective of Routing ................ 10
2.1 A System .............................................. 14
2.2 Calculation of Channel Outflow with
the System Impulse Response Function .................. 20
2.3 Development of the Analytical
Impulse Response Function of the
Linearized SaintVenant Equations .................... 27
2.4 Channel Impulse Response Function
Used to Represent Flow Hydrographs .................... 35
2.5 Analysis of Limitations of Linear Theory .............. 44
2.6 Comparison of the Inverse Gaussian and Gamma PDF
with Parameters Matched by Method of Moments
for a = 1, n = 4 and k = 2.0 ........................... 48
2.7 Comparison of the Inverse Gaussian and Gamma PDF
with Parameters Matched by Method of Moments
for a = 8, n = 8 and k = 0.5 .......................... 49
2.8 Comparison of the Inverse Gaussian and Gamma PDF
with Parameters Matched by Method of Moments
for a = 1, n = 1 and k = 0.5 .......................... 50
2.9 Development of the Linear State Space Model ........... 53
3.1 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Q0 = 10000 cfs and ain = 3 hours .......... ..63
3.2 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for 00 = 100000 cfs and ain = 3 hours ......... 64
3.3 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for = 200000 cfs and ain 3 hours......... 65
3.4 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Q0 = 10000 cfs and 0in = 6 hours .......... 66
viii
Title
A Small Inflow Transient Routed and Full SaintVenant Equations Channel for Q0 = 100000 cfs and
A Small Inflow Transient Routed and Full SaintVenant Equations Channel for Q0 = 200000 cfs and
with Linearized in a Rectangular Gin = 6 hours .........
with Linearized in a Rectangular ain = 6 hours .........
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q0 = 10000 cfs and Gin = 12 hours .........
A Small Inflow Transient Routed and Full SaintVenant Equations Channel for Q0 = 100000 cfs and
A Small Inflow Transient Routed and Full SaintVenant Equations Channel for Q = 200000 cfs and
with Linearized in a Rectangular a in = 12 hours ........
with Linearized in a Rectangular ain = 12 hours ........
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q0 = 10000 cfs and Gin = 24 hours .........
A Small Inf low Transient Routed and Full SaintVenant Equations Channel for Q0 = 100000 cfs and
A Small Inflow Transient Routed and Full SaintVenant Equations Channel for Q0 = 200000 cfs and
with Linearized in a Rectangular a in = 24 hours ........
with Linearized in a Rectangular a in = 24 hours ........
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q0 = 10000 cfs and ain = 6 hours ..........
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q0 = 100000 cfs and ain = 6 hours .........
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q0 = 200000 cfs and ain = 6 hours .........
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q 10000 cfs and a. = 12 hours .........
Chne oQ in
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q0 = 100000 cfs and in = 12 hours ........
Page
67 68 69
3.8 3.9 3.10
3.11 3.12
3.13
3.14 3.15
3.16 3.17
I
'igure Title Page
3.18 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Triangular
Channel for Q0 = 200000 cfs and ain = 12 hours ........ 81
3.19 Verification Process for the
Linear State Space Model ....... .................. .... 84
3.20 Comparison of the Linear State Space Model
and 3Parameter Gamma PDF for
a = 1, n = 4 and k = 2.0 .................................. 87
3.21 Comparison of the Linear State Space Model
and 3Parameter Gamma PDF for
a = 8, n = 8 and k = 0.5 .............................. 88
3.22 A Small Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a =in 6 hours, L = 100 miles and % = 10000 cfs ....... 92
3.23 A Small Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
( in = 3 hours, L = 100 miles and Q = 100000 cfs ...... 94
3.24 A Small Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
ain = 12 hours, L = 100 miles and % = 100000 cfs ..... 95
3.25 A Small Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a in = 24 hours, L = 100 miles and Q = 200000 cfs ..... 96
3.26 A Small Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
ain = 6 hours, L = 100 miles and Q = 100000 cfs ...... 97
3.27 A Small Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
ain = 12 hours, L = 100 miles and Q = 200000 cfs ..... 99
3.28 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a in = 3 hours and L = 100 miles ....................... 103
Figure Title Page
3.29 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a = 3 hours and L = 400 miles ....................... 104
3.30 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a in = 6 hours and L = 100 miles ....................... 105
3.31 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a i= 6 hours and L = 400 miles ....................... 107
in
3.32 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a in = 12 hours and L = 100 miles ...................... 108
3.33 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
ain = 12 hours and L = 400 miles ...................... 109
3.34 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a. = 24 hours and L = 100 miles ...................... III
in
3.35 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a in = 24 hours and L = 400 miles ...................... 112
3.36 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a. = 6 hours and L = 100 miles ....................... 113
in
3.37 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a. = 6 hours and L = 400 miles ....................... 114
in
3.38 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a. = 12 hours and L = 100 miles ...................... 115
in
x I
Figure Title Page
3.39 A Large Inflow Transient Routed with the Linear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a in = 12 hours and L = 400 miles ...................... 116
3.40 A DoublePeaked Inflow Transient Routed
with the Linear State Space Model and the
Full SaintVenant Equations in a Rectangular
Channel for L = 100 miles ......................#...... 118
3.41 A DoublePeaked Inflow Transient Routed
with the Linear State Space Model and the
Full SaintVenant Equations in a Rectangular
Channel for L = 400 miles and 0 = 100000 cfs ......... 119
3.42 A DoublePeaked Inflow Transient Routed
with the Linear State Space Model and the
Full SaintVenant Equations in a Rectangular
Channel for L = 400 miles and 1 = 0000 cfs .......... 120
3.43 A DoublePeaked Inflow Transient Routed
with the Linear State Space Model and the
Full SaintVenant Equations in a Rectangular
Channel for L = 400 miles and Q = 200000 cfs ......... 121
4.1 Approach to the Nonlinear Analysis .................... 124
4.2 Analysis of the Nonlinear Characteristics
of the Primary Model .................................. 126
4.3 Variation with Flow of the Channel
Lag for ain = 12 hours
for a Rectangular Channel ............................ 133
4.4 Variation with Flow of the Channel
Lag for Gin = 12 hours
for a Triangular Channel .............................. 134
4.5 Variation with Flow of the Channel
Standard Deviation for Gin = 12 hours
for a Rectangular Channel .............................. 135
4.6 Variation with Flow of the Channel
Standard Deviation for Gin = 12 hours
for a Triangular Channel .............................. 136
4.7 Variation with Flow of the Channel
Coefficient of Variation for Gin = 12 hours
for a Rectangular Channel ............................. 138
'igure Title Page
4.8 Variation with Flow of the Channel
Coefficient of Variation for in = 12 hours
for a Triangular Channel .............................. 139
4.9 Approach to the Development of the
Nonlinear State Space Model ........................... 143
4.10 A Cascade of Nonlinear Reservoirs
Governed by Sm = K0 .................................. 144
4.11 Development of the Nonlinear
State Space Model ....... .1........... ........ 147
4.12 A General Strategy to Emulate
Any Primary Model ..... ........ ......... ........... 159
5.1 Verification Process for the
Nonlinear State Space Model.......................... 161
5.2 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
ain = 3 hours and L = 100 miles....................... 165
5.3 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
0 in = 3 hours and L = 400 miles ....................... 168
5.4 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
Gin = 6 hours and L = 100 miles.................... 170
5.5 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
Cin = 6 hours and L = 400 miles ....................... 171
5.6 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
a. = 12 hours and L = 100 miles ...................... 172
in
5.7 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
U in = 12 hours and L = 400 miles ...................... 173
44
Figure Title Page
5.8 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
ain = 24 hours and L = 100 miles ...................... 175
5.9 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Rectangular Channel for
ain = 24 hours and L = 400 miles ...................... 176
5.10 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a in = 6 hours and L = 100 miles ....................... 177
5.11 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a in = 6 hours and L = 400 miles ....................... 178
5.12 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a in = 12 hours and L = 100 miles ...................... 180
5.13 A Large Inflow Transient Routed with the Nonlinear
State Space Model and the Full SaintVenant
Equations in a Triangular Channel for
a in = 12 hours and L = 400 miles ...................... 181
5.14 A DoublePeaked Inflow Transient Routed
with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular
Channel for L = 100 miles ....... o ................... 183
5.15 A DoublePeaked Inflow Transient Routed
with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular
Channel for L = 400 miles and 0 = 100000 cfs ......... 184
10
5.16 A DoublePeaked Inflow Transient Routed
with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular
Channel for L = 400 miles and 0 = 10000 cfs .......... 185
5.17 A DoublePeaked Inflow Transient Routed
with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular
Channel for L = 400 miles and 0 = 200000 cfs ......... 186
How to Emulate Any Primary Model ...................... 194
 1 11
figure Title Page
6.2 A Piecewise Constant Versus Flow
Relationship for m ...... .............................. 198
7.1 Lag Versus Mean Discharge Relationship,
South Creek Experimental Catchment,
University of New South Wales ......................... 204
7.2 A Comparison of the Nonlinear State Space and
MuskingumCunge Models with the Full SaintVenant
Equations at L = 400 miles in a Rectangular Channel... 207
Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy A NEW PERSPECTIVE OF
NONLINEAR ROUTING IN HYDROLOGY
By
George Francis Smith
December 1983
Chairman: Wayne C. Huber
Major Department: Environmental Engineering Sciences
Existing routing models for hydrologic systems have either 1) a complex, nonlinear structure not amenable to second order (i.e., error) analysis, requiring considerable computational resources for solution; or 2) a simple, linear structure that is strictly valid over only a small flow range about which the model is linearized.
The new perspective for routing surface runoff or channel flow uses any appropriate model, called the primary model, to simulate the observed physics for the catchment or channel of interest, and then emulates the mathematical behavior of the primary model with a computationally efficient nonlinear state space model that has structural properties amenable to second order analysis. The overall strategy matches impulse response properties of the state space model to those of the primary model over a wide range of flow levels. In general, impulse response properties vary with flow level.
xvi
Two important response properties are time lag and attenuation.
Nonlinearities in time lag occur because water flows faster as it flows deeper. Nonlinearities in attenuation occur because deep waves disperse more rapidly than shallow waves. These two properties can be expressed as linear or piecewise linear functions of discharge in logspace, thus relating them directly to the parameters of a nonlinear state space model.
The limitations of this approach to routing were examined by
comparing results from the nonlinear state space model with a numerical solution of the SaintVenant equations. No significant differences were observed. The state space model required two to three orders of magnitude less CPU time than numerical solution of the SaintVenant equations.
A generalized, computationally efficient state space routing model that can emulate the nonlinear behavior of any other routing model for downstream waves has been developed. Further studies are required to determine the full usefulness of this modelling approach.
Xvi i
CHAPTER I
A NEW PERSPECTIVE OF NONLINEAR ROUTING
Statement of the Problem
Existing technology for routing unsteady flow in open channels and for routing surface runoff in a catchmient has not proved adequate for some important hydrologic applications. Current routing methods are either inefficient when numerous simulations are required, as in flood frequency analysis; inadequate when the forecast requires an estimate of uncertainty in past or future flow conditions, as in water supply forecasting; or unsatisfactory when the information content of current data is important, as in realtime river forecasting. The above applications require
1. an objective estimate of uncertainty,
2. repetitive computations, and/or
3. consideration of second order (i.e., error) properties of
models and data as well as the first order (i.e., mean)
trajectory in time and space of the flow conditions;
as well as requiring models or computational algorithms to account for significant nonlinearities of the channel or catchment response.
For example, the National Weather Service (NWS), the Bureau of Reclamation, and the Federal Emergency Management Authority typically study the thousands of dams throughout the country by simulating various dam failure scenarios. Numerical solution of the full SaintVenant equations can impose a significant computational burden in such
contingency dam break simulations for which repetitive calculations are required. Simpler models, while computationally efficient, do not properly account for the nonlinear behavior important to these analyses.
The NWS studies of various possible dam failures are designed to provide timely and accurate forecasts of expected flow levels to communities downstream of the dam in the event of a dam failure. When the threat of a dam failure is imminent, the NWS Forecast Office can provide worst case (i.e., instantaneous failure) and alternative projected flow conditions to the media and emergency authorities. Several computer simulations involving flow routing downstream from each dam studied are required to produce that information.
For daily, realtime river forecasting, as performed at NWS River Forecast Centers, timely processing of all current observations is critical. A rapid mechanism is needed to update the conditions projected by the computer model to account for the most recently observed data. A state space model of the full SaintVenant equations, which would be amenable to updating with modern estimation theory, imposes a large computational burden. Simpler routing models are not hydraulically adequate over the range of inflow transients encountered in flood forecasting. Although some simple models could be formulated for use with updating algorithms, their forecasting capability is diminished~by the limited flow range over which they properly model the nonlinear nature of water movement.
Flood frequency studies, such as those performed during the hydraulic design of channel structures, often require computer simulations for many events over a wide range of flow conditions. Again, numerical solution of the full SaintVenant equations is costly
and the simpler, linear models do not properly simulate for large variations of inflow.
Statement of the New Perspectv
The new perspective of routing in hydrology is to let any model, called the primary model, simulate the physics of a channel or catchment, and to emulate the behavior of the primary model with a model structured to meet the applications described above. This process is depicted in Figure 1.1. To account for the modelling of second order properties (i.e., variance and covariance measurements of uncertainty) and to simplify computation, a state space form has been chosen as the model to emulate primary model behavior. The state space model structure is
X =F X + G 'tU(1) where Xt = the system state vector,
.Et= the state transition matrix,
2GX,t = the input coefficient matrix, and
U+l= the vector of system inputs.
The structure of a state space model is possibly the only one that meets the previously described objective of simplicity and has the potential for simulating second order properties. The state space model can be used with modern estimation theory in a Kalman filtering algorithm to propagate and update both mean and variance properties of the channel or catchment. The primary model simulates the physics of the hydrologic system while the state space model, because of its desired computational and structural properties, emulates the mathematical behavior of the primary model.
PRIMARY MODEL OF THE OBSERVED PHYSICS
NONLINEAR STATE SPACE
MODEL
EMULATE THE MATHEMATICAL
BEHAVIOR OF THE PRIMARY MOCEL
Figure 1.1
An Overview of the New Perspective of Nonlinear Routing in Hydrology
The essence of this new perspective of routing is its twostep approach to the problems encountered with existing routing technology. First, the physics of the system are modelled with a primary model, then the mathematics of the primary model are emulated by a state space model. The state space model is a model of the primary model.
Development of the New Perspective
Overview
Several steps are required to apply this new perspective to a routing problem of interest. The steps are
1. select a primary model,
2. calibrate the primary model to the hydrologic system,
3. determine the appropriate structure for the state space
model,
4. calibrate the state space model to the primary model, and
5. emulate the primary model behavior with the state space
model.
Selection of a Primary Model
The primary model simulates the physics of a channel or
catchment. Because, in general, the forces governing the motion of water interact in a complicated manner, a complex model typically will be most appropriate to describe the physical properties of the hydrologic system. However, the primary model does not have to be complex. Any simple model can be chosen to simulate the system physics based on data or computational limitations, or modeller preference. A
summary of some available channel routing models is presented in Appendix A.
Calibration of the Primary Model to the Hydrologic System
The effort required to calibrate the primary model will vary
greatly depending upon the complexity of the primary model chosen. If a physically detailed hydraulic model is selected, channel crosssection data must be obtained and reduced to the form required by the model. This requires a complicated process of adjusting friction factors, determining active versus inactive flow area, etc. For simpler models, the calibration process usually requires determining a few parameters based on the comparison of simulated and observed data for some past inflows. Parameter selection for simple models can often be automated, given inflow and outflow records. Even for some complex hydraulic models, final parameter adjustment can be performed automatically (Fread and Smith 1978).
Structure of the State Space Model
The mathematical properties of the primary model can be
represented in numerous ways, given the general structure of equation
1.1. The system states (i.e., the elements of the vector X) can represent many possible quantities. The states of the surrogate model are determined by the important relationships in the primary model. State space models of hydrologic systems can be constructed with physically based states such as the quantity of water in a reservoir (SzollosiNagy 1981) or with states which are purely mathematical
artifacts that reproduce observed system behavior (Goldstein and Larimore 1980). In fact, the full SaintVenant equations can be expressed in state space form. At each time step, values of flow and crosssectional area at numerous points along the channel represent the state of the system. The state vector, X, could be constructed from these flow and area values. However, since the number of points defined along a channel typically is large, the state vector would be correspondingly large. For example, some NWS applications of the Dynamic Wave Operational Program (DWOPER), which numerically solves the SaintVenant equations (Fread and Smith 1978) have 25 and 45 points defined along the channels, respectively, for the ColumbiaWillamette River system and the OhioMississippi River junction. The state vectors for these systems would consist of 50 and 90 elements, respectively, with state transition matrices of 2500 and 8100 elements. Obvious dimensionality problems arise in the computer solution of equation 1.1 for systems of these sizes.
The structure of the nonlinear state space model developed is presented in detail in Chapter IV. The approach taken uses the nonstationary linear structure of the state space model as a surrogate for the stationary nonlinear behavior of the primary model.
Calibration of the State Space Model to the Primary Model
Problems encountered in certain routing applications can he solved by using the state space model to emulate the mathematical behavior of the primary model. Since the state space model is calibrated to the mathematical properties of the primary model, many approaches are possible.
The approach taken here draws on the general, powerful theory of linear systems. As described in Chapter II, the nearly linear behavior of the primary model for a small flow range is emulated with an equivalent linear system. This is accomplished by matching the impulse responses of the primary and state space models for small flow ranges about various reference discharges.
The inherent nonlinearities that cause the behavior of the state space model to agree with or to depart from the primary model behavior are variations with discharge of the lag and dispersive properties of the impulse response function. These properties vary substantially over large changes of inflow. A simple functional relationship can be derived that specifies the variation of F and G with discharge to emulate properly the lag and dispersive properties. Given such a functional relationship, the state space model can emulate the full nonlinear and dynamic behavior of the primary model over all flow ranges.
The State Space Model as Surrogate for the Primary Model
Simple linear models have been used in place of complex routing models for many years. They replace the complex model, but do not properly represent some of the important properties of the complex model. With the new perspective, the primary model is not replaced. instead its mathematical structure is emulated by a state space model with computational and structural properties appropriate for the solution of problems encountered in many routing applications. The state space model acts as a surrogate for the primary model.
The state space model typically has much lower dimensionality than the primary model and imposes a much smaller computational burden. Because the state space model emulates the full nonlinear behavior of the primary model, it can be used in applications where the primary model is not appropriate for error analysis because of the computational burden or improper structure. This approach presents a model mathematically equivalent to the primary model, that can be used on smaller computers than would be possible with the primary model itself.
An Application of the New Perspective
To explore the new perspective's applicability to routing in hydrology, the general steps described above were followed for a specific primary model. As shown in Figure 1.2, the application selected was the full, nonlinear SaintVenant equations as the exemplar primary model to simulate downstream flow in a prismatic channel. This primary model and physical system were chosen because
1. the SaintVenant equations are generally recognized as the
most complete representation of unsteady flow phenomena,
2. an analytical solution for the impulse response of the
SaintVenant equations can be derived for downstream wave
movement in prismatic channels, and
3. numerical solution of the SaintVenant equations is
wellstudied, and many excellent solution algorithms are
available for comparison with the surrogate model (Liggett
and Cunge 1975).
Analysis of the ability of a state space model to emulate the behavior of the full, nonlinear SaintVenant equations should lend
INEW PERSPECTIVE
OF ROUTING I
SELECT SRINTVENANT EQUATIONS
AS EXEMPLAR PRIMARY MODEL
N
DEVELOP LINEAR DEVELOP NONLINERR
STATE SPACE MODEL STATE SPACE MODEL
VERIFICATION OF NONLINEAR STATE
SPACE MODEL
7f
7/
GENERAL PROCEDURE TO
EMULATE ANY PRIMARY MODEL
Figure 1.2 An Overview of the Development
of the New Perspective of Routing
LINEAR RNRLYSIS OF SAINTVENANT
EQURT IONS
NONLINEAR ANALYSIS OF SRINTVENANT EQURTATIONS
VERIFICATION OF LINEAR STATE SPRCE MODEL
I
D
N
insight into the applicability of this approach for other primary models. Because, in general, no analytical solution for the impulse response of the primary model is available, the case chosen allows comparison of the true impulse response with results obtained numerically.
Organization of the Dissertation
The organization of this dissertation is schematically presented in Figure 1.2. The new perspective of routing and the reasons for selecting the SaintVenant equations as the exemplar primary model have been presented. A linear analysis of the SaintVenant equations and development of a linear state space model follow in Chapter II. In Chapter III, results from the linear state space model are compared with the full SaintVenant equations for small inflow transients about a reference flow level. These results indicate that the linear state space model is a very adequate surrogate for the SaintVenant equations over small flow ranges. Further nonlinear studies are motivated by comparison of the linear state space model with the primary model for large inflow hydrographs. Nonlinear behavior of the primary model causes the comparison of results with the linear state space model to deteriorate.
A nonlinear analysis of the SaintVenant equations primary model and the development of a nonlinear state space model are conducted in Chapter IV. Chapter V is devoted to the verification of results for the nonlinear state space model versus the primary model for large inflow hydrographs. Results from the primary and surrogate models for a complex inflow hydrograph were comparable. However, the nonlinear state
space model used two to three orders of magnitude less CPU time than the numerical solution of the SaintVenant equations.
A general procedure to emulate any primary model with the
nonlinear state space model is presented in Chapter VI. By following the steps here described, the nonlinear state space model can be used as a surrogate for any channel routing or surface runoff model.
The relationship of the current work to previous studies of routing in hydrology is the subject of Chapter VII. A review of pertinent work by other authors, some data to support the power function relationship between flow level and time lag, and a comparison of the nonlinear state space and MuskingumCunge models with the full SaintVenant equations are the major topics of this chapter.
Conclusions on the new perspective of routing in general, and on its application to the exemplar primary model in particular, are stated in Chapter VIII. Recommendations for future work, including the necessary tasks for the application of estimation theory with the nonlinear state space model and possible alternative approaches for representing the nonlinear model structure, form the final portion of this last chapter.
CHAPTER II
INSIGHTS INTO NONLINEAR ROUTING FROM LINEAR SYSTEMS THEORY Linear Systems Theory
Introduction
Linear systems theory provides many tools for the analysis of
linear phenomena. Linear systems theory, strictly applied, would have few uses because, in general, the only truly linear systems are mathematical abstractions. However, the spirit of linear systems analysis lends itself to unraveling those aspects of nonlinear systems that behave in a nearly linear manner. The phenomenon of routing in hydrology is, in general, a very nonlinear process. Nevertheless, important insights into hydrologic system behavior may be gained by applying some of the tools of linear systems theory.
The concepts of linear systems theory have been applied to many and varied systems and are well known. Those concepts germane to this work are summarized below.
Systems can be thought of abstractly as the mechanisms that interrelate two objects (Dooge 1973, pp. 34). The operational definition used in this work is that a system is the process that transforms input signals into output signals. For the hydrologic processes analyzed here, the complex relationships among the forces governing the motion of particles of water are conceptualized in a system as shown in Figure 2.1. The value of this view of a system lies
INPUT SIGNALS
x1(t) x2(t)
SIGNALS
.(t)
Figure 2.1 A System
SYSTEM
OUTPUT
in the analytical tools it provides since the mathematical properties of the system determine which tools are applicable.
Properties of a System
The properties of consequence for hydrologic systems are whether the system is
a. determinsitic or probabilistic,
b. distributed input/distributed output
C. causal,
d. stationary, and
e. linear.
Typically the response of a timevarying, nonstationary, distributed, probabilistic, nonlinear hydrologic system is simulated by a time invariant, stationary, lumped, deterministic, linear model. And we wonder why there are discrepancies between observations and simulated outputs!
A system is determinsitic if the precise state of the system (the values of X in equation 1.1) can be foretold. If a random component is introduced into the variation of the system states, the system becomes probabilistic.
The inputs and outputs of most hydrologic systems do not interact at a single point. As such, there should properly be some spatial description of the system processes. The complex, dynamic, nonlinear models of channel routing of surface runoff typically account for some of the spatial variation of parameters over the channel or catchment. The parameters of simple linear models usually apply over the entire channel or catchment area and, therefore, such models operate As if all
input/output interactions occurred at a single point (i.e., they are lumped models).
In general, a system can have several inputs and outputs as shown in Figure 2.1. The hydrologic systems analyzed in this work are assumed to have a single input and a single output.
A causal system is one in which outputs cannot precede inputs.
All physical systems are causal. This restriction limits the range over which integral equations describing the system are applied.
The stationarity of a system is defined by whether or not the
input/output relationship changes with time. If a system is stationary, identical input produces identical output.
Linearity is the property that opens the doors to the storehouse of mathematical tools. If system inputs are defined as xj, system outputs as yip and the transformation of x into y by the system is symbolized as > then given
x >y and x
a system is said to be linear if
a I xI + a 2 x2 > a I yI + a 2 y2 (2.1)
where a 1 and a2 are constants (Thomas 1969, p. 140).
The analysis of nonlinear systems is much more difficult than that of linear systems because the principle of superposition applies to linear systems (Thomas 1969, p. 139). Superposition, which follows directly from the definition of linearity, says that if inputs are added or scaled, outputs are similarly added or scaled. Therefore, the system behavior can be analyzed independently of the inputs. This is a point of great significance: if the system output for an instantaneous unit
input can be found, the system behavior for any input can be obtained by considering the actual input as a series of scaled unit inputs at various instants of time (Liu and Liu 1975, pp. 135146). Given the system output for an instantaneous unit input (called the impulse response), the output for a stationary linear system can be obtained for any input with the convolution integral
t
y(t) = f h(tT)'x(T) dT (2.2)
where x('), y(.) are the system inputs and outputs respectively, and
h(') is the system impulse response (Johnson and Johnson 1975,
pp. 137140).
Stationary nonlinear behavior can sometimes be simulated by a nonstationary linear model. This is, in fact, exactly the strategy behind the nonlinear state space model developed in Chapter IV. The theoretical basis is most clearly seen by recognizing that for a linear model, the input/output relationship can be represented by matrix multiplication as
[IT,,r," "Ij" h h h hn
i1 2 3 m 11 12 13 in
0 h h . .
22 23
0 0 h .
33 (2.3)
0 0 0 . h
mn
[0,02, 0 O 0n]
1 2
where Ii = the input at time i, an element of row vector I, Oj = the output at time j, an element of row vector 0, and
hi,j = an element of the input/output relationship matrix H.
18
Note that I and 0 are time vectors, in contrast to the state vector X of equation 1.1. Matrix H transforms I into 0. What we want at any time, t, is to generate one entry, Ott in vector 0. To do this we need the appropriate weights to apply to all inputs Ii up to and including It. These weights will fill one column of matrix H. However, it is physically impossible to observe the weighting function of a system. An approach which circumvents that problem is to generate the model impulse response. This is equivalent to solving for 0 in equation 2.3 when I is all zeroes except for a single entry of 1 in, say, the ith location. Call this vector 6 The output, 0, generated will be the ith row of matrix H. The input, 6,, selects the ith row of H and presents it as the output, 0. For a linear nonstationary system the rows of H will, in general, be different from each other (i.e., each row, it will represent the system impulse response at time i). The assumption implicit when using a single system impulse response in the convolution integral is that the system is linear and stationary. En this case, the ith row of H will be equal to the (il)th row of H shifted one position to the right. Therefore, equation 2.3 can be rewritten for a stationary linear system as
I~ h h h h h
 1 2 3 4 ~ m
o h h h .
1 2 3
o 0 h h .
1 2
o 0 0 h 0 (2.4)
o 0 0 h 1h2
0 0 0 0 h
The desired weighting functions (i.e., the columns of H) are seen in equation 2.4 to be the timereversed images of the system impulse response for a stationary linear system. Figure 2.2 shows how the convolution integral (equation 2.2) exploits this to compute the system output.
Nonlinear system behavior can be modelled by letting the columns of H in equation 2.3 vary each time step. The elements on any diagonal will not be equal, as is the case for the linear stationary system represented by equation 2.4. Therefore, the system impulse response (the rows of H) will change each time step, thus permitting a nonstationary linear model to account for stationary nonlinear behavior. The state space model developed in Chapter IV emulates the nonlinear system behavior by varying the weighting function based on the level of flow.
A General Linear Differential Equation
Linear timeinvariant equations can be used to approximate the behavior of many physical systems. This is a good approximation when system characteristics change very slowly relative to variations of the inputs (Ogata 1967, p. 307). Channel and catchment systems can be well represented by linear timeinvariant equations.
A general differential equation for a nonstationary linear system can be written as
a (t) + a (t)"dnly(t) +
n n ni ni
dt dt
(2.5)
+at)dy(t)
+ aW) dt + a0(t).y(t) = x(t)
20
TIMEREVERSED IMAGE
OF THE SYSTEM
IMPULSE RESPONSE
TIME INTO THE PAST k]
EL.
A I
,'t (PRESENT)
THE CHANNEL OUTFLOW AT ANY TIME IS A WEIGHTED AVERAGE OF THE PAST AND
PRESENT INFLOWS. THE TIMEREVERSED IMAGE OF THE SYSTEM IMPULSE RESPONSE
REPRESENTS THE WEIGHTING FUNCTION.
(t) hk). I tdc o 0 .
N, ( PRESENT )
Figure 2.2 Calculation of Channel Outflow with
the System Impulse Response Function
where x(t) = the system input,
y(t) = the system output,
a0(t), al(t), . ., anl(t), an(t) are coefficients, and
the initial conditions are specified for (t n1dnlY(to)
Y(t) . and dn
0 dt "1dt
The system is linear because there are no powers or products of y. Equation 2.5 can be rewritten in matrix form as
d _(Y(t)) A(t)*Y(t) = B(t)*X(t) dt
Y(t) =
A(t) =
y(t) dy(t)
dnly(t)
dtn1
0
0
0 0
Sa0(t) al(t)
a
I
B(t) = I
n(t)
0
0
0 1
a n1(t a (t) ni
a~ (t)
, and
(2.6)
where
0
x(t) =
0
x(t)
If the system is stationary the coefficients A(t) and B(t) of the underlying differential equation 2.5 are not functions of time. In this case equation 2.6 can be rewritten as
(Y(t)) A.Y(t) = BX(t) (2.7)
dt
which can be solved using the integrating factor etA to yield
Y(t) = e(tt)Y(t) + f e(t)AB.X(r ) dt (2.8)
to
(Ogata 1967, p. 315).
The first term on the right hand side of equation 2.8 is the zero input response of the system. With no input the system would move toward equilibrium starting from this point when t = t0. The second term is the zero state response, which describes the motion of the system driven by input X(T) for t 0< T < t.
Transformation of Moments by a Stationary Linear System
One way to describe a system's input and output signals is to study their moments. The following section, on the use of transform methods to study moments of stationary linear systems, is based on work by Dooge (1973, pp. 113115).
The Rth moment of a function about the time origin is
Uj(f) = f f(t)'tR dt (2.9)
and the Rth moment about the center of mass is
U(f) f(t)(tU')R dt (2.10)
The output signal of a stationary linear system can be expressed with the convolution integral (equation 2.2) in terms of the input signal and the system impulse response. The relationships between moments of the input, output and impulse response are most easily seen by transforming equation 2.2 with either the Fourier or Laplace transform. Since moments about the center of mass of the input and output signals are of interest, we will use the bilateral Laplace transform defined by
FB(s) f f(t).est dt (2.11)
where f(t) is a function in the time domain.
Differentiating equation 2.11 R times and setting s = 0 produces the following expression for the Rth moment about the origin in terms of the Laplace transform (Dooge 1973, p. 114).
R
U (f) = (1)R 2,FB(S)j (.=0
Rds RB s=0
The convolution integral can be expressed in the transform domain by taking the bilateral Laplace transform of equation 2.2 to yield
YB (s) HB(s)YXB(s) (2.13)
where YB(S), HB(s) and XB(s) are the Laplace transforms of y(t), h(t)
and x(t), respectively (Dooge 1973, p. 114).
The Rth moment of the output about the origin follows from equation 2.12
as
R.
U,(y) = (1)R'AR YB(s)]s=O (2.14)
dsR
Substituting equation 2.13 for YB(S) gives
UR(y) = ( )R' H (s)XB(s)]s=O (2.15)
ds
Solving equation 2.15 with Leibnitz's rule for continued differentiation of a product, and substituting for the moments with equation 2.12 yield
R
UR(y) = ()*URk(h)*Uk(x) (2.16)
k=0
where (R) indicates a combination of R things taken k at a time, or
k
( ) R!
k k!*(Rk)!'
U' (h) = the (Rk)th moment about the origin of the system
Rk
impulse response, and U'(x) = the kth moment about the origin of the system input
k
signal (Dooge 1973, p. 114). If, instead of about the origin, the moments are taken about the centers of mass of the input, output and system impulse response functions, the following relationship holds (Dooge 1973, p. 115)
R
UR(Y) = ()*URk(h)*Uk(x) (2.17)
k=0
Only the first moment about the origin and the second moment about the center of mass are needed for the ensuing work. For R = 1 equation 2.16 reduces to
(2.18)
U'(y) = U'(h) + U'(x)
1 1 1
Since U &)=0, for R = 2 equation 2.17 becomes
U 2(y) U 2(h) + U 2(x) (2.19)
Therefore, for the first moment about the origin (the time lag) and the second moment about the center of mass (the variance), the moment of the system output is the sum of the moments of the system impulse response and the system input (Dooge 1973, p. 115).
Summary
The impulse response is the output generated by an instantaneous unit input to a system. Because the impulse response function completely describes the behavior of a stationary linear system, we would like to determine the impulse response of the channel or catchment physical systems. The best description of the real world that we have is the primary model. The goal of linear analysis of the system is to determine the impulse response of the primary model at various reference flow levels. The changing properties of the impulse response function with flow level lend insight into the nonlinear behavior of the physical system.
An operational way to determine the impulse response of the
primary model is to pulse the model with a unit input and observe the output. That is, in general, the approach that will be taken to find a system's impulse response.. An analytical solution can be found for the impulse response of the SaintVenant equations of unsteady flow for downstream waves in a prismatic channel. To derive the impulse response, the SaintVenant equations must be linearized about a reference discharge. This process is described In the next section.
Linearization of the SaintVenant Equations for a General Prismatic Channel
To determine analytically the impulse response for routing
unsteady flow, the SaintVenant equations must be linearized in terms of small perturbations about a reference flow condition. The steps presented in this section are outlined in Figure 2.3. Dooge (1973, pp. 245253) presented the special case for linearizing the SaintVenant equations for a rectangular channel and Chezy friction formula. The ensuing derivation follows Schaake (1980) who extended Dooge's work for a general prismatic channel and resistance formula.
The SaintVenant equations for unsteady flow can be expressed as
3A a0 ay a(Av)
+ B. t+ = q (2.20)
and
g + ax+ a + S S =q V) (2.21)
g X '9 x f 0 gA (Ux
where A = cross sectional area,
Q = discharge,
t = time,
x = distance along the channel,
v = average velocity,
y depth of flow,
B = surface width of flow,
q = lateral inflow per unit width,
ux = the x component of the velocity of the lateral inflow,
g = the acceleration of gravity,
Sf = the friction slope, and
LINEARIZE ABOUT FLOW LEVEL,
REPRESENT
CHANNEL GEOMETRY AS
B BONA'
SRINTVENANT
EQUATIONS
Fi~q IT
I
IMPULSE RESPONSE Is INVERSE GAUSSIAN PDF
Figure 2.3
Development of the Analytical Impulse Response Function of the Linearized SaintVenant Equations
LINEAR,
SECOND ORDER HYPERBOLIC
EQUATON
REPLACE HIGH ORDER
TIME DERIVATIVES
WITH SPACE
DERIVATIVES USING
KINEMATIC WAVE
CELERITY
LINEAR,
SECOND ORDER PARABOL I C EQUATION
so = the channel bottom slope (Henderson 1966, p. 287).
Equation 2.20 is the continuity equation for a channel. The effects of the forces acting on a control volume within the channel are modelled by equation 2.21, the momentum equation. Since the continuity equation is already linear, only the momentum equation requires linearization.
Equation 2.21 is first rewritten in terms of the dependent
variables Q and A. To this end, the friction term and channel geometry are expressed as
V = k.W 7P. (2.22)
and
B =B A r (2.23)
0
where B is the top width,
B0 is the top width for flow QO2 and
k, p and r are described below.
Equation 2.22 is a general relationship for friction in open channel flow which can be related to the wellknown Manning or Chezy friction formula. With the approximation
R A A (2.24)
h P B
where Rh =the hydraulic radius and
P =the wetted perimeter.
Equation 2.22 becomes Chezy's friction formula with
k C
and (2.25)
p =1/2
or Manning's friction formula with
k = 1.49/n
and (2.26)
p = 2/3
where C = the Chezy flow resistance factor and
n = the Manning flow resistence factor (Daily and Harleman
1966, pp. 297298).
Equation 2.23 can approximate most regular cross sectional shapes. For a rectangular channel r = 0; while for a triangular channel r = 1/2. However, a trapezoidal channel is not well modelled by equation 2.23 for all flow ranges. For shallow flow a trapezoidal channel behaves similarly to a rectangular channel; while for deep flow it can be approximated by a triangular channel. There is a transition flow range where equation 2.23 does not represent the change in top width with area for a trapezoidal channel. Substituting 2.23 into 2.22 and rearranging gives B2Q
Sf 0 (2.27)
k 2 A2[p(lr)+l]
Let
m = p(lr)+l (2.29)
and
a 0 (2.29)
BP 0
where m is analogous to the dimensionless kinematic wave parameter defined by Eagleson (1970, p. 250). Substitution of 2.28 and 2.29 into 2.27 gives
S Q2 S (2.30)
f a2 A2m 0
Substituting equations 2.20 and 2.30 into the momentum equation gives
3 A 23Q'"Q 2Q
(g A3_B*Q 2) 3x + 2*B QA + BA 2. =
(2.31)
gBA3+rS gBQ2 3+r2m gB*A *S ___ *
0 2
where henceforth the terms involving q will be ignored because we are studying prismatic channels with no lateral inflows. Finally, differentiating 2.20 with respect to x, differentiating 2.31 with respect to t while holding the lefthand coefficients constant, and
equating the common @2A term, results in
g*A v2 2Q 2Q 2Q
(gA v2) 2*v*
B 3x2 3x*8Tt 2
(2.32)
2mg*S +2 S 0" x + v at
Linearizing equation 2.32 about Q0 A0, v0 and B0 forms the linear, second order, hyperbolic equation
(g A v)'.2Q 22v 2Q
B 0 a 2 0 FxP t 8t2
(2.33)
2*g*S
2 m g S *"3 Q + v 0 t
0
where A0, v0 and B0 are the cross sectional area, velocity and top
width, respectively, for the reference flow, % (Schaake 1980).
The impulse response of the linearized form of equation 2.33,
presented in Appendix A, is a complex function containing a Dirac delta and a modified Bessel function (Harley 1967, p. 18). A simpler, approximate solution can be found with an assumption from kinematic wave theory. The 3x.3t and 3t2 terms of equation 2.33 can be replaced by a ax2 term using the relationship
a c" (2.34)
where
dx
c (the wave celerity).
This substitution is routinely made in the derivation of diffusion routing equations (Koussis 1976). Here its use results in the simplified linear approximation to the complete equations
3Q DQ Q 320
 + Co.  2 B 0 (142) (2.35)
at ax 2 S0x2
where
2 = F2. [1m(2m)] (2.36)
0
in which
Q /A
F = 0 0 (the Froude number at Q ) (2.37)
0 Vg.A /B 0
0 0
and
c = m'v (the kinematic wave velocity) (Schaake 1980). (2.38)
0 0
The diffusion form of equation 2.35 appears often in channel and pollutant transport modelling (Fischer 1967, Carslaw and Jaeger 1959, Bansal 1971). Koussis (1976) followed a derivation similar to the one
presented above, but was interested in a simpler numerical solution to the SaintVenant equations, not an analytical one. The interest in this section is in an analytical solution to equation 2.35. This will provide the theoretical basis for comparison of results from a numerical solution of the full SaintVenant equations, and the linear and nonlinear state space models developed later in this chapter and in Chapter IV, respectively.
The influence of all terms in the SaintVenant equations has been retained in the parabolic form of equation 2.35. The impulse response function of equation 2.35 is
(xc t)2
h(x,t) = x exp{ 4. Dt } (2.39)
4. .D. t3
where
0
D = 2B 142) (2.40)
2B S
0 0
(Harley 1966, Schaake 1980). The impulse response, h(x,t), defines the flow for all positive x and t. The area under h(x,t) is unity, (i.e., f h(x,t) dt = 1) to preserve continuity, and the units are
0
time1. Equation 2.39, the solution for an instantaneous input at x = 0, t = 0 in equation 2.33, appears in many scientific disciplines as a solution to the advectiondiffusion equation (Carslaw and Jaeger 1959, p. 50). Since the diffusivity, D, is constant for any reference flow, Q0, Fickian diffusion is represented by equation 2.39 (Slade 1968, pp. 8081).
Insights from Linear Systems Theory
Introduction
The analytically derived impulse response for the SaintVenant equations is strictly correct only for small variations about the reference flow. In the spirit of linear approximations for nonlinear systems, the impulse response can be used to approximate the channel response for all flow ranges, as is shown in Chapter III.
Equation 2.39 completely describes the response of the general
prismatic channel for flows close to Q,. This analytic solution can be manipulated to gain insights into the behavior of the SaintVenant equations, from which 2.39 was derived.
The Gaussian form of the solution of equation 2.39 for fixed t is well known (Crank 1956, pp. 911). However, Schaake (1980) recognized that for fixed x the form of equation 2.39 is an inverse Gaussian probability density function (pdf). Given this insight, a wealth of information can be obtained from previous studies of this statistical distribution (Johnson and Kotz 1970, pp. 137153).
The Inverse Gaussian PDF as an Analytical Solution of the Unsteady Flow Equations
Since 1871, when Barre De SaintVenant proposed the equations of unsteady flow in an open channel, people have sought an analytical solution for equations 2.20 and 2.21. No analytical solution to the full SaintVenant equations has been found. An analytical expression for the outflow hydrograph would be ideal for studying the behavior of a channel. However, the ideal has not been attained because, in addition
to the lack of an analytical solution of the SaintVenant equations, there is no analytical expression for a general inflow hydrograph.
There is, however, an analytical solution to the equations of
unsteady flow for a special class of inflow hydrographs in a prismatic channel. Because equation 2.39, the equation for an inverse Gaussian pdf, is the response for an impulse input to a prismatic channel, an inverse Gaussian pdf can be observed at any point along the channel. Parameters of the inverse Gaussian pdf are determined by the channel properties and the distance, x, from the point at which the impulse is input. As shown in Figure 2.4, an inverse Gaussian pdf with x = x0+ L is the analytical expression for the outflow hydrograph when an inflow hydrograph is specified by an inverse Gaussian pdf with x = x0. Channel properties determine the parameters c and D of equation 2.39. This means that for the class of inflow hydrographs defined by inverse Gaussian pdfs, the outflow hydrographs from a prismatic channel are also inverse Gaussian pdfs (Schaake 1980)! This is of great significance because the parameters of the outflow inverse Gaussian pdf are related simply to the parameters of the inflow inverse Gaussian pdf and the channel properties.
Moments of the Inverse Gaussian PDF
The time delay and dispersion characteristics of 2.39 are
expressed by its first two moments (Johnson and Kotz 1970, pp. 137140). The first moment about the origin (i.e., the mean) describes the time lag properties as
k t = (2.41)
UPSTREAM IMPULSE
RT x 0O
IMPULSE RESPONSE
fT x xo
I
IMPULSE RESPONSE
RT x xo+ L
LFICTICIOUS CHRNNELI
LENGTH x.
XIXO xxo
REAL CHANNEL
LENGTH L
xxo+ L
Figure 2.4 Channel Impulse Response Function
Used to Represent Flow Hydrographs
x0
The second moment about the mean (i.e., the variance) describes the dispersion properties as
Q "x
k = o2 Q 0 (142) (2.42)
2 B .S .c3
0 0 0
The coefficient of variation, a dimensionless measure of the impulse response, can be obtained from 2.41 and 2.42 as 1/2
I 0 I
c = B 0 1,(I_2)I (2.43)
v B 0 S 0 c x
The form of the coefficient of variation, equation 2.43, indicates that the kinematic wave assumption is most nearly correct for shallow
Q
flow (smallB *Q ) on steep slopes (large SO). Woolhiser and Liggett
B *c
0 0
(1967) define a dimensionless parameter to determine when the kinematic wave assumption is appropriate
L S
kWL 0 0 0 (2.44)
WL H F2
0 0
where L0 characteristic length,
So = channel bottom slope,
H= depth of flow, and
F2 = Froude number.
0
If kWL > 10, the SaintVenant equations are approximated well by a kinematic wave solution (Eagleson 1970, p. 336). These conditions are generally met by overland flow, with its shallow depth of flow and steep slopes. For river channels with steep slopes, the So term of the momentum equation dominates, making the kinematic wave assumption a good
approximation of the SaintVenant equations. For rivers with flat slopes, the SO and terms are dominant in equation 2.21 (Henderson 1966, p. 364). The diffusion analogy model, discussed in Appendix A, uses only these two terms for channel routing.
The Inverse Gaussian PDF as an Inflow Hydrograph
The interrelationships between peak flow (Qp), celerity (cO) and timing properties (as measured by the variance, a2 ) of the inflow in
hydrograph can be expressed for a given total flow volume (I) (Schaake 1980). With these relationships several properties of hydrographs input to routing models can be controlled. Any three of Qp' co, 2n or V can be specified arbitrarily. Since equation 2.39 is the solution for a pulse input at x = 0 and t = 0, we will define a point downstream at x = x0 to observe the flow for t > 0. With this notation and h(*,*) defined by equation 2.39, the total flow volume, inflow hydrograph, and peak discharge are, respectively,
= f Q(x 0,t) dt = f V.h(x 0,t) dt = Vf h(x 0,t) dt (2.45)
0 0 0
Q(x ,t) = V.h(x ,t) (2.46)
0 0
Q = V'h(x ,t (x )) (2.47)
where h(x,t) is a maximum when
3D c *x 2 2
t (x) = 1 + 1 (2.48)
0
38
By choosing any three of V, Qp, Gin and co, a hydrograph with
certain desired properties can be selected. Since equation 2.39 defines the response to an impulse input, an inverse Gaussian shape hydrograph can be observed at any point along the channel. The channel can be divided into two sections as shown in Figure 2.4. The upstream, ficticious section has an impulse input and produces an inverse Gaussian pdf outflow hydrograph which becomes the inflow hydrograph for the downstream, real channel of interest. When an impulse is input to the upstream section at x = 0, an inverse Gaussian pdf appears as input to the downstream section at x = xO. The variance of the flow transient as x = x0 (jin) can be selected by proper choice of x. Letting a2 and x of equation 2.42 equal a2 and x0, respectively, gives in 0 epciey ie
Qx
a2 Q 0 0 (1_02) (2.49)
in B S 'c3
0 0 0
Substituting equation 2.40 into equation 2.49 simplifies the expression for the inflow variance to
2 D* x
C2 0 (2.50)
in c2
0
The parameters D and cO are specified by the channel geometry and slope. The value of x can be arbitrarily selected to give the desired inflow variance. Equation 2.50 can be rearranged to solve for x0 as
c3.a2
x = __0 in (2.51)
0 2"D
39
The channel section where 0 < x < x is a ficticious reach designed solely to generate an inflow hydrograph with desired properties for the real channel being modelled. This real channel is defined by the x0 < < x + L channel section where L is the length of the real channel. The outflow hydrograph is also an inverse Gaussian shape with variance
2"D(x + L)
out= 3 (2.52)
out c3
0
Equations 2.18 and 2.19 can be rewritten to solve for the mean and variance of the impulse response
U'(h) = U'(y) U'(x) (2.53)
and
U (h) = U (y) U (x) (2.54)
2 2 2
With these relationships and equations 2.41 and 2.50, the mean of the impulse response function for the linearized SaintVenant equations can be expressed as
x+L x
U'(h) = 0 0 (2.55)
1 c c
c0 c0
or
U'(h) = lag L (2.56)
1 channel
and the variance of the impulse response function is
2. D'(x + L) 2D'x U (h) = a2 G2 0 0 (2.57)
2 out in c3 c3
0 0
2.DL (2.*58)
U (h) = C7 oL( .8
2 channel c3
0
Equations 2.56 and 2.58 express the lag and variance, respectively, which are added to an input hydrograph as it passes through a prismatic channel specified by parameters D, c0 and L. The channel impulse response will vary depending on the reference flow level chosen, since D and c0 are functions of Q0. Note, however, that the channel impulse response is independent of x0. This means that once QO is chosen, the impulse response of the system is not dependent on the input signal.
The development of a linear state space model which follows later in this chapter uses the properties of the linearized SaintVenant equations described above. The remainder of this section follows Schaake (1980) and focuses on several properties of equation 2.39 that lend additional insight to routing in hydrology.
Additional Properties of the Inverse Gaussian PDF
Attenuation of the peak inflow, Qp, to the real section of the channel can also be determined from an analysis of equation 2.39. Define the attenuation of the peak inflow at a point distance L below the upstream location, x0, as
Y.h(x + Lt (x + Q)
A(L) 0 (2.59)
V.h(x ,tp (x 0)
where the numerator is the peak outflow rate at location x0+ L and the denominator is Qp. Cancelling the common factor V produces htx + L,tp(x +t)
A(L) = 0 p 0 (2.60)
h(x 0, tp (x 0))
Schaake (1980) notes the interesting result that although A, the inflow hydrograph attenuation, is independent of V, the downstream peak outflow rate is directly proportional to V.
The channel lag of the center of mass (i.e., the mean) of the
inflow hydrograph is given by equation 2.56. The peak lag time can be derived from equation 2.39 as
tp= tp(x0+ L) tp(X0) (2.61)
or 1/2 1/2
tp + 0 + 0 (2.62)
=X _3'D_, D 39A Linear State Space Approach
Introduction
The inverse Gaussian pdf is a linear approximation of the full nonlinear SaintVenant equations. For small variations in flow the inverse Gaussian pdf produces results nearly identical with the full equations (Chapter III). Unfortunately, derivation of a state space model for equation 2.39 is cumbersome. A state space model with parameters related to the parameters of the inverse Gaussian pdf is derived in this section.
The inverse Gaussian pdf is defined by equation 2.39 in terms of the parameters co) the wave velocity, and D, a diffusion constant (,Johnson and Kotz 1970, p. 137). This form of the inverse Gaussian pdf follows naturally from the physically based parameters of the
SaintVenant equations. Equations 2.38 and 2.40 define co and D, respectively, in terms of the reference discharge and cross sectional properties. However, in the ensuing analysis, the system impulse response will be determined from its moment properties which will be deduced from the changing moments of the inflow and outflow, using equations 2.56 and 2.58. To facilitate the determination of the inverse Gaussian pdf from its moment properties, an alternative form of equation
2.39 can be defined in terms of two different parameters: p, the time mean of the distribution, and X, an inverse measure of dispersion (Johnson and Kotz 1970, p. 138). Equation 2.39 can be rewritten in terms of P and X as
h(x,t) : flG(t1',x) =X/(2 I.t3),exp{X'(tV)2/(2"2"t)} (2.63)
The parameters of equations 2.39 and 2.63 are related by
x (2.64)
C
0
S2(2.65)
D
The mean (PIG), variance (OG), and skewness (YIG) of the inverse Gaussian pdf are simple functions of the parameters p and X (Johnson and Kotz 1970, pp. 139140).
IG (2.66)
02 = P3/) (2.67)
IG
Y = 23 6)
(2.68)
These simple relationships facilitate the determination of state space model parameters. This is demonstrated in the following sections.
Criteria for Matching Impulse Response Functions
The value of a linear system analysis of the unsteady flow
equations is that is allows the development of a simple model that has an impulse response function similar to the impulse response of the full SaintVenant equations. For small perturbations about a reference flow, Q0, equation 2.63 defines a system impulse response which is the same as the impulse response for the full SaintVenant equations. In Chapter III, we will see that equation 2.63 is an excellent approximation to the SaintVenant equations for small (10 percent of Q0) inflow transients. The linear approximation deteriorates for large (2000 percent of Q) inflow hydrographs. The results are also presented in Chapter III to help show the limitations of the linear theory. The strategy employed to determine the limitations of the linear theory is depicted in Figure
2.5. The primary thesis of this work is that the fundamental nonlinearities present when routing large inflow hydrographs with the SaintVenant equations can be accounted for in a model that matches the impulse response of the SaintVenant equations at all reference flow levels, Q0"
There are a number of possible criteria which could be used to compare impulse response functions for the simple and complete models. In general, functional relationships for the impulse responses cannot be obtained, so the criteria for matching impulse response functions must be able to compare the responses of the systems to impulse inputs. Some of the methods available for comparison of impulse responses are
SRINTVENRNT EQUAT IONS
ESTIMATION OF
IMPULSE RESPONSE FUNCTION PARAMETERS
HHm
LINERRIZED SRINTVENANT EQUATIONS
ANALYTICAL IMPULSE RESPONSE
FUNCTION
COMPARISON OF IMPULSE RESPONSE
FUNCTION
UNDERSTANDING OF THE LIMITRTICNS CF THE LINEAR THEORY
Figure 2.5
Analysis of Limitations of Linear Theory
TRANSIENT NUMERICAL
INFLOW SOLUTION HYDROGRAPHS (DWOPER)
OUTFLOW
HYDROGRAPH (NUMERICAL)
1. least squares regression,
2. maximum likelihood analysis,
3. spectral analysis, and
4. matching moments.
Direct derivation of parameters for the least squares criterion may not be easy. For the simplest of models, a functional relationship may be found between the parameters and the criterion to be minimized. Differentiation of that function for each parameter and solution of the set of the simultaneous equations generated is almost always difficult because of nonlinearities (Dooge 1970, pp. 170171). Statistical packages are available to perform least squares regression and to determine parameter values but, for any except the simplest models, those packages can be costly to use.
Maximum likelihood parameter identification suffers from the same weaknesses. It may be difficult to determine the likelihood function for nonlinear models and the cost of using statistical packages to identify parameters involves significant cost for complex models.
Spectral analysis, which involves transformations into the
frequency domain, comparison of selected properties, and transformation back into the time domain, in principle, can be carried out to determine parameter values. The state space and full systems can be equated by matching pole and zero locations of the transfer functions in the transform domain. Since most hydrologic situations are not described functionally, the transformations to and from the spectral domain must be carried out numerically. Numerical inversion of the transform back to the time domain is almost always difficult (Dooge 1973, p. 31). The autocovariance properties of the impulse response function previously
have been used to find the parameters of a state space model for a unit hydrograph (Goldstein and Larimore 1980). This approach was not pursued, but may prove useful for determining state space model parameters when the impulse response is not a singlepeaked function.
The mean and variance of the impulse response functions for the
linear and full systems can be compared to determine parameter values by the method of moments. If the first n moments of two impulse responses are identical, the two systems will produce identical output for polynomial input of degree n or less (Dooge 1970, pp. 170171). The method of moments is used for subsequent analysis of the state space model because it is simple and allows the convenient parameterization of the model in terms of the moments of the impulse response function.
Computing the Parameters of a Gamma PDF
A cascade of linear storage reservoirs is a wellstudied routing model (Zoch 1934, Nash 1959, Dooge 1973, p. 176). This model is amenable to state space formulation (SzollosiNagy 1981) and can be functionally represented by a 3parameter gamma pdf (Johnson and Kotz 1970, p. 166)
G(t) = k'((ta).k)nl'exp{(ta)k} (2.69)
r (n)
where a = the pure delay,
n = the number of reservoirs,
k = the time delay in each reservoir, and
r() = the gamma function (Abramowitz and Stegun 1964, p. 255).
As formulated in equation 2.69, n can take on any positive
value. The interpretation of n as the number of reservoirs restricts it
to positive integers, implying that r(n) could be replaced by (nl)! (the factorial function) in 2.69. Since it is useful to work with continuous parameters, the positive integer restriction on n will be relaxed temporarily. The concept of integer values for n will be reinstated for implementation of the state space model.
Equations for the mean (PG), variance (02), and skewness (YG) of the 3parameter gamma pdf are
PG = a + n/k (2.70)
2 = n/k2 (2.71)
G
Y = 2// (2.72)
(Johnson and Kotz 1970, pp. 167168). Parameters of the inverse Gaussian and gamma pdfs can be related by equating their means, variances and skewnesses as defined by equations 2.66 through 2.68 and 2.70 through 2.72. When the gamma pdf parameters are found in terms of the inverse Gaussian parameters, the surprisingly simple results are
1
a = U (2.73)
4 X
n = 4 X (2.74)
k = 2 *X (2.75)
T *U2
As demonstrated in Figures 2.6 and 2.7, an inverse Gaussian pdf can be approximated well with a gamma pdf using equations 2.73 through
2.75. The relationships based on matching moments begin to deteriorate as n approaches 1, because the gamma pdf becomes simply a shifted exponential (Johnson and Kotz 1970, p. 166). Results for this case are shown in Figure 2.8.
3PARAMETER GAMMA POF CASHED LINE WITH SYMBOLS
a 1.0 n 4.0 k 2.0
INVERSE GAUSSIAN POF
SOLID LINE p 3.00 x 27.00
0.0 1.0 2.0 3.0 4.0
TIME
Figure 2.6
5.0 6.0 7.0 8.0
Comparison of the Inverse Gaussian and Gamma PDFs with Parameters Matched by Method of Moments fora=1,n=4andk=2.0
3PARAMETER GAMMA PDOF DASHED LINE WITH SYMBOLS
a 8.0 n 8.0 k 0.5
INVERSE GAUSSIAN POF
SOLID LINE p 24.00 x 432.00
50.0
Figure 2.7
TIME
Comparison of the Inverse Gaussian and Gamma PDFs with Parameters Matched by Method of Moments for a = 8, n = 8 and k = 0.5
C
0
o
0
a,
3PARAMETER GAMMA POF DASHED LINE WITH SYMBOLS
a 1.0 n 1.0 k 0.5
INVERSE GAUSSIRN POF SOLID LINE p 3.00 x 6.75
0.0
Figure 2.8
2.0 4.0 6.0 8.0 10.0
TIME
Comparison of the Inverse Gaussian and Gamma PDFs with Parameters Matched by Method of Moments for a = 1, n= 1 and k = 0.5
f
I'
I I'
I'
I
II
0I
/:
Qi
c :
12.0
Parameter values for a cascade of linear reservoirs model can be computed from values of P and X with equations 2.73 through 2.75. Values of P and X can be determined from channel properties with equations 2.64 and 2.65 if x, c0 and D are known, or from the mean and variance of a channel response function with
= IG (2.76)
W3
IG (2.77)
U2
IG
(Johnson and Kotz 1970, pp. 139140). Equations 2.76 and 2.77 can be substituted into equations 2.73 through 2.75 to provide relationships for the parameters of the gamma pdf in terms of the mean and variance of the system impulse response function (i.e., the inverse Gaussian pdf) as
1
a = IG (2.78)
n=4IG (2.79)
IG
k =.* I (2.80)
G a2
IG
In practice the mean (II) and variance (aG would be estimated from the moments of the primary model input and output hydrographs. To be strictly correct in this case, the inverse Gaussian mean and variance in equations 2.78 through 2.80 should be replaced by p2 and o2 IG IG'
respectively, where the indicates a value estimated from data.
The goal of this approach is to emulate a calibrated primary routing model (such as the numerical solution of the SaintVenant equations) with a linear state space model. The general approach is to
1. simulate with the primary model using an inverse Gaussian
pdf as the input hydrograph,
2. compute the mean and variance of the input and output
hydrographs,
3. estimate the mean (YIG) and variance (a2G of the system
impulse response with equation 2.56 and 2.58, and
4. compute a, n and k from equations 2.78 through 2.80.
When the values of a, n and k computed in step 4 above are
substituted into equation 2.69, G(t) defines the impulse response of the desired linear system. A state space model with an impulse response identical to equation 2.69 is derived in the next section.
A Cascade of Linear Reservoirs State Space Model
The state space model developed in this section is based
conceptually on a cascade of linear reservoirs. Figure 2.9 provides an overview of the steps taken to develop the linear state space model. The relationship among inflow, outflow and storage for a reservoir is described by the lumped continuity equation
dS = 0 (2.91)
dt
where S = the reservoir storage,
I = the reservoir inflow, and
0 = the reservoir outflow.
A reservoir is linear when the relationship between storage and outflow is linear as
S = .0 (2.82)
where K = a constant.
GENERAL STATE SPACE
EQUATION
4i FX +G tUi
I
I
RESERVOIR CONTINUITY EQUATION
AS/At 1 0
Figure 2.9 Development of the Linear
State Space Model
LINEAR RESERVOIR EQUATION S C.0
CHAN
RESPO
AT
FL
NEL IMPULSE NSE FUNCTION REFERENCE OW LEVEL
EQUATIONS FOR LINEAR RESERVOIRS
IN A CASCADE
EXPRESS F AND GIN TERMS OF a, n AND k
COMPUTE a,n,k AT REFERENCE FLOW LEVEL
STATE SPACE MODEL WITH
F AND G AS FUNCTIONS
OF CHANNEL IMPULSE RESPONSE
I
I/
The state space model will operate with discrete time steps of At when implemented, so rewriting equation 2.81 for discrete time yields
AS (2.83)
At
where I = (1P)*I + P*1 1 2
0 = (1p)*0 + p*0
1 2
Subscripts 1 and 2 indicate the beginning and end of a time step, respectively, and p can take on values from 0 to 1; T and O are weighted average values of the inflows and outflows, respectively, at the beginning and end of a time step. Substituting for I and 0 in equation 2.83 gives
AS
AS = (1p)* I + p*I (p) p.02 (2.84)
At 1 2 1 2
Rewriting AS in terms of storage at the beginning and end of a time step and substituting for S with equation 2.82 produce
K0 K*0
2 1 = (1p)*I + p.I (1p).O0 p.O0 (2.85)
At 1 2 1 2
2 2 1 2
(2.86)
K*0 (1p)0 *At
1 1
or, solving for 02
A< (1p)*At]
0 = At (1p)I + p*I ) + ( p0t) (2.87)
2 ( + PAt) 1 2 + pAt) 1
Equation 2.87 applies for any reservoir in a cascade. The form of equation 2.87 varies slightly for the first reservoir because the
inflow, It, is externally specified for the first reservoir. For all subsequent reservoirs, the inflow is equal to the outflow from the previous reservoir, as
Iti = 0t'iI ; i > 1 (2.88)
where the first subscript indicates the time step and the second subscript indicates the relative position of the reservoir in the cascade.
Rewriting equation 2.87 with this double subscript notation for the first reservoir yields
(K (1P)At)
t+1 (()I + p t+1) + (2.89)
where It is the inflow to the first reservoir at time t and
= K + p.At.
The governing equation for all subsequent reservoirs is
Ot+1,i = ((uPo'tiI + P'Ot+i1I +
(2.90)
(K (IP)A t)
0ti
where i > 1.
Recall the general form of a state space model from equation 1.1
X =F "X +G U
,+ X t t ,'t t+1
Equations 2.89 and 2.90 must be organized into the matrix form indicated by equation 1.1. For a linear system the state transition and input coefficient matrices are constant, so equation 1.1 can be simplified to
X t+l = FX + GU (2.91)
Defining the state vector for a cascade of n linear reservoirs as
X =
0t 1
Ot,l
t,2
0
t,n
and the system input vector as
U
t+ 1
(1P)* It + t+1
0 0
(2.93)
F and G can be determined by substituting equations 2.89, 2.90, 2.92 and
2.93 into equation 2.91. Defining
S= (a (1p)*At)
At
S= (1p) + p.
(2.94)
(2.95) (2.96)
the F and G matrices for a cascade of n identical linear reservoirs are
S0 0 0
T 0 *
F = p*T2*G T*D T 0 (2.97)
n2 nI 2
p .*T Q ~ pT .Q T
(2.92)
and
T
G:T2 (2.98)
nl n
Parameters k (equation 2.80) and K (equation 2.94) are related by
K = i/k (2.99)
The parameter a (equation 2.69) is accounted for in the state space model by delaying the time at which the outflow from the nth reservoir appears by a units of time.
Summary
A linear state space model that emulates the behavior of the full SaintVenant equations about a reference flow level has been developed. Parameters of the state space model are determined by matching the mean and variance of the system impulse responses of the primary and surrogate models. Verification of the cascade of linear reservoirs state space model by comparing the results with equation 2.39 (the impulse response of the linearized SaintVenant equations) and a study of the limitations of the linear state space model for large inflow hydrographs are the major topics of Chapter III. The work of Chapter IV, in which the important nonlinearities are analyzed and a nonlinear state space model is developed, is motivated by the problems arising when large inflow hydrographs are routed with the linear state space model.
CHAPTER III
LIMITATIONS OF THE LINEAR THEORY Introduction
The linear theory developed in Chapter II for prismatic channels is applicable strictly only for small perturbations around a reference flow level. In this chapter, outflows predicted by the linearized SaintVenant equations and the linear state space model are compared with numerical solutions of the full SaintVenant equations for both small (10 percent of baseflow) and large (2000 percent of baseflow) inflow hydrographs.
The model used as the numerical solution of the SaintVenant equations is the National Weather Service Dynamic Wave Operational (DWOPER) Model (Fread 1978). This dynamic wave routing model is based on an implicit finite difference solution of the complete onedimensional SaintVenant equations of unsteady flow.
a(A+A )
a at q = 0 (3.1)
ax B t+Q 3(Q2/A) + .A(h
3 + (Q2/A) gA( h+ Sf) q'v + Wf.B = 0 (3.2)
Bt 8x fxx x
where Q = discharge,
A = cross sectional area,
A0 = offchannel storage area, q = lateral inflow or outflow,
x = distance along the channel,
t = time,
g = the acceleration of gravity, h = the water surface elevation,
vX = velocity of lateral inflow in the xdirection,
Wf = the wind term,
B = channel top width, and
Sf = the friction slope defined as
s n .o (3.3)
2.2.A2. R4/T
in which n = the Manning roughness coefficient,
R = the hydraulic radius, and
1" = the absolute value function (Fread 1978).
Offchannel storage, lateral inflows and wind effects are ignored in this work.
The solution technique in the DWOPER model is implicit; therefore, the time step size can be selected based on accuracy rather than numerical stability. This makes the DWOPER model very efficient in the use of computer time. Simulation of actual river reaches typically requires from 10 to 30 seconds of CPU time on an IBM 360/195 computer (Fread 1978).
Verification of the Linear SaintVenant Equations
The impulse response function of the linearized SaintVenant equations was derived in Chapter Il as equation 2.39.
h(x,t) =,exp{ Lt2
V 44n Dt TDt}
where
Q
D T 0 s 142)
0 0
and
C = mev 0 0
In this section, results obtained by simulating with the DWOPER model and by solving equation 2.39 are presented for 10 percent of baseflow transients on baseflows of 10000, 100000 and 200000 cfs. Inverse Gaussian shape hydrographs with Gin = 3, 6, 12 and 24 hours were routed through a rectangular channel (width = 100 feet) with the DWOPER model. Transients for the same three baseflows, but for Gin = 6 and 12 hours were routed through a triangular channel with side slopes of 1:1. The bottom slope for both channels was 1 foot per mile.
Comparisons between the numerical solution of the full
SaintVenant equations and the analytical solution of equation 2.39 were made at 100 and 400 miles from the upstream end of the channels. The values of co and D in equation 2.39 were found from the channel properties at each baseflow level. Table 3.1 shows the combinations of baseflow, Gin and channel geometry that produce the results for small inflow transients.
It should be emphasized that we are comparing results from the numerical solution of a set of quasilinear partial differential equations with results from the analytical solution of the inverse Gaussian pdf. The area under the curve defined by equation 2.39 will be the same for all values of x. The area under a hydrograph is the volume of water in the transient. Because there are no gains or losses of water in the prismatic channels studied, the volume must be equal for all x. This is guaranteed by equation 2.39.
Organization of Figures 3.1 through 3.18, Comparing the Linearized and Full SaintVenant Equations
Rectangular Channel
Gin (hours)
6 1 12
10 Figure 3.1 Reference
Flow 100 Figure 3.2
(1000 cfs)
200 Figure 3.3
Figure 3.4 Figure 3.5 Figure 3.6
Figure 3.7 Figure 3.8 Figure 3.9
Figure 3.10 Figure 3.11 Figure 3.12
Reference
Flow
(1000 cfs)
10
100
Triangular Channel
in (hours)
6 1 12
Figure 3.16 Figure 3.17 Figure 3.18
Table 3.1
3 
i I
Figure 3.13 Figure 3.14 Figure 3.15
I
I I
Figures 3.1, 3.2 and 3.3 show results for Gin = 3 hours for baseflows of 10000, 100000 and 200000 cfs, respectively. Outflows simulated with the DWOPER model and the solution of the inverse Gaussian pdf agree closely with respect to the timing of the hydrographs for 100 and 400 miles from the channel inlet, but differ in the magnitude of the flows. In these cases the discrepancies are caused by problems with the numerical solution of the SaintVenant equations. For an abrupt wave, such as the Gin = 3 hour inflow hydrograph, the distance and time step sizes for the numerical solution must be small. For the results presented in Figures 3.1 through 3.3, the distance step was 2.5 miles and the time step was 0.25 hours. With so many distance and time steps over 400 miles and 100 hours, numerical roundoff caused the DWOPER model to not maintain continuity. Therefore the flows predicted by the DWOPER model differ from those found with equation 2.39. However, the timing properties of the very abrupt inflow hydrograph are excellent.
Slightly less severe inflow hydrographs, Gin = 6 hours, are
depicted in Figures 3.4, 3.5 and 3.6. Results from the inverse Gaussian pdf and the DWOPER model compare very well, especially for the higher flows. The numerical roundoff problems of the DWOPER model are much less evident for this inflow hydrograph.
For the longer duration inflow, with ain = 12 hours, presented in Figures 3.7, 3.8 and 3.9, the numerical solution of the full SaintVenant equations is matched well by the inverse Gaussian pdf. For these results, as well as those shown in Figures 3.10, 3.11 and 3.12 for Gin = 24 hours, the analytical and numerical solutions are almOst indistinguishable. For these broader hydrographs the distance and time step sizes can be increased so that the numerical roundoff is no longer a problem.
10.8
10.64
I.
10.4i
10.2i
:.s
* ...
LEGEND
INFLOW HYDROGRAPH
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES 1 IMPULSE RESPONSE 400 MILES
100
0
0
0
%.0
Cd
C9
cr
C)
o
a
M 75 19 r 'r"N 'Fv .
0 10 20 30 40 50 60 70 80 90
TIME (HOURS)
Figure 3. 1 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 10000 cfs and a, = 3 hours
Im M
LEGEND
INFLOW HYDROGRAPH
........................... ..........................
o DWOPER MODEL 100 MILES
o DHOPER MODEL 400 MILES
 IPL REPONSE 100 MILES
* IMPULSE RESPONSE 400 MILES a IMPULSE RESPONSE 400 MILES
5 10 15 20 25 30 35 40 45 50 55 TIME (HOURS)
Figure 3. 2 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 100000 cfs and a, = 3 hours
110108
106
104
102
* .
* .
* .
* *
* .
* .
* .
* *
~~1
.601W
LEGEND
INFLOW HYDROGRRPH....
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
* IMPULSE RESPONSE 400 MILES
5 10
40 45
I 25 :
TIME (HOURS)
Figure 3. 3 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 200000 cfs and a, = 3 hours
220
216
212208 
* is'
I I'
204 200
.
LEGEND
INFLOW HYDROGRAPH
..... ao.~o.oo ..... ....... .. o. ol.o~ o
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
 NS  4d
* IMPULSE RESPONSE 100 MILES
* IMPULSE RESPONSE 400 MILES
20 30 40 50 60 70 80 90 100 110 120 TIME (HOURS)
3. 4 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 10000 cfs and a, = 6 hours
10.8
10.6
10.4
10.2 1
I
j
0 10
Figure
130
130
ja
LEGEND
......... INFLOW ...HYDROGRAPH ..........................
o DWOPER MODEL 100 MILES o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES a IMPULSE RESPONSE 400 MILES
10 20 30
10 20 30
40 50
TIME (HOURS)
Figure 3. 5 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 100000 cfs and a, = 6 hours
110]
1081
106
104 
102
I
. I m ~mf
m
LEGEND
......... INFLOW HYDROGRAPH .............
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
* IMPULSE RESPONSE 400 MILES
2M
1 10 20
0 10 20
Figure 3. 6
* I
* S.
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Qo = 200000 cfs and a, = 6 hours
40 50
TIME (HOURS)
216
212 208204
200
LEGEND
INFLOW HYDROGRAPH
o DWOPER MODEL 100 MILES
a DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
10.2
10
Figure 3. 7
60 80 100 120 140 160 180
TIME (HOURS)
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Qo = 10000 cfs and a, = 12 hours
11
10.8
10.6 10.4
Il0a
I I Ir I I I
10 20 30 40 50 60 70 80 90
TIME (HOURS)
1 12 110 110 120 130
 100 MILES
 400 MILES
 100 MILES
.FL. LEGEND
................. .........8........
! o DWOPER MODEL
! 1 DWOER MDELDWOPER MODEL
S\ IMPULSE RESPONSE
Ss IMPULSE RESPONSE
*
! i
S.*..
 '***
140 150 160
Figure 3. 8
A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Qo = 100000 cfs and a, = 12 hours
 400 MILES
10oe8
106I
104 
102100
r
416
216
216
GEND
ORAPH
I". .....""T i'[
 100 MILES
 400 MILES PONSE 100 MILES PONSE 400 MILES
20 30 40 50 60 70 80 0 100 110 120 130 140 TIME (HOURS)
Figure 3.9 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 200000 cfs and a = 12 hours
*.. LE
/: INFLOW HYDR
S \ o DWOPER MODE
DHOPER OCE
* IMPULSE RES
! IMPULSE RES
II I........ 
2121
204
10A
10
150
 m m I a 
LEGEND
INFLOW HYDROGRRPH
'............ ....................
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
1110.8
250 275 300 325 350
TIME
375 400 (HOURS)
425 450 475 500 525 550
Figure 3.10 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 10000 cfs and a, = 24 hours
a IMPULSE RESPONSE 400 MILES
10.4
10.2
200
I
LEGEND
INFLOW HYDROGRAPH
o...............+ee. ~o ....................... l~l+ o ....o+.. l...
o DWOPER MODEL 100 MILES
o DHOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
* IMPULSE RESPONSE 400 MILES
Irn
 I I I I
125 150 175 200 225
TIME
i
250
(HOURS)
27 0O I I I
275 300 325 350 375 400
Figure 3.11
A Small and Full Channel
Inflow Transient Routed with Linearized SaintVenant Equations in a Rectangular for Qo = 100000 cfs and a' = 24 hours
110
108
106104 
102
100
220 216
o
0 212
0
ce CS
~2084
T
54
200
75
I.L LEGENDD
......... INFLOW HYDROGRAPH...........
/o DWOPER NODEL 100 MILES
[] DWOPER MODEL 400 MILES IMPULSE RESPONSE 400 MILES
! IMPULSE RESPONSE 400 MILES
I R*
: S
* \
* S
! S
;Q
! i
* 5
\.
3 125 150 175 200
Figure
225
250I
w~in m .
75
275
300
325
TIME (HOURS)
3.12 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Rectangular
Channel for Qo = 200000 cfs and a, = 24 hours
To demonstrate the ability of the inverse Gaussian pdf to match the DWOPER model for a general prismatic channel shape, defined by equation 2.23, results for ain = 6 hours in a triangular channel are presented in Figures 3.13, 3.14 and 3.15. The delay in the arrival of the peak flow and the increased dispersion caused by the triangular channel shape can be seen by comparing, for example, Figure 3.4 with Figure 3.13. Simply by changing co and D of equation 2.39 to account for the new channel shape, the inverse Gaussian pdf remains an excellent approximation to the DWOPER model when the inflow transient is small.
Figures 3.16, 3.17 and 3.18, for ain = 12 hours in a triangular
channel, can be compared with Figures 3.7 through 3.9. Again the timing and dispersive properties of the channel are identified well by the inverse Gaussian pdf.
Since the impulse response function of a linear system is a
function of the system not of the inputs, co and D have the same values for a given baseflow for all rectangular channel results. Similarly, for the triangular channel, co and D are functions only of baseflow, not of ain It is helpful in the analysis of, for example, Figures 3.1,
3.4, 3.7 and 3.10 (each of which is for Q = 10000 cfs in a rectangular channel), to realize that only x0, as defined in Figure 2.4, is changing. The sets of results for the same baseflow level for either the rectangular or triangular channel shape can be interpreted as an impulse response at everincreasing distances from x = 0, the location of the upstream impulse.
Figures 3.1 through 3.18 demonstrate that, for a wide range of flow levels and hydrograph durations in prismatic channels defined by equation 2.23, the impulse response function of the linearized
119w j
C .:
10.8
C,
1 "
0 i a 0.6 :
a:
C) i :
* 1
1o.2 :
In
LEGEND
INFLOW HYDROGRAPH
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILE
* IMPULSE RESPONSE 400 MILES
(N
60 70 80 TIME (HOURS)
Figure 3.13 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Triangular
Channel for Qo = 10000 cfs and o = 6 hours
ISO
LEGEND
INFLOW HYDROGRAPH
o DWOPER MODEL 100 MILES
o OHOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES m IMPULSE RESPONSE 400 MILES
40 50 60
TIME (HOURS
110108
Figure 3.14 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Triangular
Channel for Qo = 100000 cfs and a, = 6 hours
106104 
vEr
p I 0 9 F I 9
70 80 90 100 110 120
5)
220
216 1
0
C) 212 :
0
C)
r.I.
CD
"I I
2W4
204'
200
0
LEGEND
INFLOW HYDROGRPH
o WOPER MODEL 100 MILES............
a DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES IMPULSE RESPONSE 400 MiLs
10 20 30 40 50
TIME
60
(HOURS)
80 90 too00
Figure 3.15 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Triangular
Channel for Qo = 200000 cfs and o, = 6 hours
110
A
LEGEND
......... INFLOW HYDROGRAPH ..........................
o DWOPER MODEL 100 MILES
a DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
11
10.8 
10.2
10
1 20 0 20
S40
40
Figure 3.16
" '
60
A Small and Full Channel
80 100 120 140 160 180u
TIME (HOURS)
Inflow Transient Routed with Linearized SaintVenant Equations in a Triangular for Qo = 10000 cfs and a,, = 12 hours
10.6 4
10.4
200
200
I
m IMPULSE RESPONSE 400 MILES
LEGEND
INFLOW HYOROGRRPH
......................ffc .................
o DWOPER MODEL 100 MILES
o DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
I I I
40 60 80o
TIME (HOURS)
110
100
Figure 3.17
A Small and Full Channel
Inflow Transient Routed with Linearized SaintVenant Equations in a Triangular for Qo = 100000 cfs and a, = 12 hours
108
106 104 
102
160
m IMPULSE RESPONSE 400 MILES
LEGEND
......... INFLOW HYDROGR PH...............
o DWOPER MODEL 100 MILES
a DWOPER MODEL 400 MILES
* IMPULSE RESPONSE 100 MILES
* IMPUI.SE RESPONSE 400 MILES
216 212 2081 204 200
0
110 120
130 1;0 150 160
Figure 3.18 A Small Inflow Transient Routed with Linearized
and Full SaintVenant Equations in a Triangular
Channel for Qo = 200000 cfs and a. = 12 hours
2201
:
* S
* S
10 20 30 40 50 60 70 80 90 TIME (HOURS)
SaintVenant equations (i.e., the inverse Gaussian pdf) produces results almost identical to those from the DWOPER model for small inflow transients. Since results from the inverse Gaussian pdf are found by analytically solving equation 2.39, the CPU time required is small (i.e., less than a second). The CPU times for the numerical solution of the SaintVenant equations on a PRIME 750 computer are presented in Table 3.2. The DWOPER model used from I to 3 orders of magnitude more CPU time than did the impulse response function to predict the flows shown in Figures 3.1 through 3.18. In the next section the impulse response of the linear state space model is compared with the inverse Gaussian pdf.
Verification of the Linear State Space Model for Small Inflow Transients
Comparison with the 3Parameter Gamma PDF
The three parameters of the linear state space model specified by equations 2.91 through 2.98 are a, the pure delay; n, the number of reservoirs; and k, the time delay in each reservoir. The parameters of the 3parameter gamma pdf (equation 2.69) can be interpreted identically to the parameters of the state space model. In this section, results from the linear state space model and the 3parameter gamma pdf are compared when the same values of a, n and k are used as parameters for both the state space model and the pdf.
Figure 3.19 depicts the strategy used to verify the linear state
space model. The parameters co and D of the inverse Gaussian pdf can be determined for a specified prismatic channel with known slope and shape
CPU Times for the Numerical Solution of the SaintVenant Equations on a PRIME 750 Computer for Small Inflow Transients
Rectangular Channel
Yin (hours)
Reference
Flow
(1000 cfs)
10 100
200
3 I
6 1
12 1
I F290.15 sec 156.65 sec 49.39 sec 20.90 sec 215.23 sec 137.00 sec 40.50 sec 15.82 sec 182.50 sec 127.08 sec 39.05 sec 15.75 sec
Triangular Channel
Gin (hours)
6 1 12
Reference
Flow
(1000 cfs)
10 100
200
144.56 sec 141.76 sec 128.87 sec
Table 3.2
39.16 sec 30.74 sec 30.66 sec

Full Text 
PAGE 1
A NEW PERSPECTIVE OF NONLINEAR ROUTING IN HYDROLOGY BY GEORGE FRANCIS SMITH A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1983
PAGE 2
ACKNOWLEDGMENTS This study was made possible by the National Weather Service Fulltime University Assignment Program, which funded the author's studies at the University of Florida for the 19781979 academic year. The Director of the NWS Hydrologic Research Laboratory, a position formerly held by Dr. Eugene Peck and more recently by Dr. Michael Hudlow, has provided initial encouragement and continuing support of the work presented here. Dr. John Schaake, Chief of the NWS Hydrologic Services Division, has guided the author's research during the past several years. His support and always carefully considered advice on both the technical and personal aspects of an effort of this magnitude are gratefully appreciated. Dr. Wayne Huber, Professor of Environmental Engineering Sciences, University of Florida, has performed masterfully the difficult task of committee chairman for absentee doctoral research. His thoughtful review of this material and concern for both the overall technical aspects of the work and the details of presentation have substantially improved the final dissertation. The comments and exchange of ideas provided by Dr. Konstantine Georgakakos, NOAANRC Research Associate, have been invaluable to the author's continuing research. The author's colleagues in the NWS Hydrologic Research Laboratory, especially Mr. Gerald Day, have provided encouragement and stimulating discussions throughout this endeavor. Their continued expressions of concern are deserving of acknowledgment. ii
PAGE 3
The careful reviews and editorial comments of Ms. Ellen Stockdale and Melanie Quinn have improved immensely the quality of this report. The computer simulations and figures were generated on the NWS Office of Hydrology PRIME 750 computer. The author thanks Mr. Gary Kemp, of the NWS Hydrologic Services Division, for his assistance in maintaining the computer support needed for this work. For their care and diligence with the final computer plotting of the figures, the author is especially grateful to Ms. Cynthia Brown, of the NWS Hydrologic Services Division, and to Mr. Stephen Ambrose and all the technical support staff of the NWS Hydrologic Research Laboratory. The help of Ms. Janice Lewis, of the NWS Hydrologic Research Laboratory, in the use of the NWS Dynamic Wave Operational Program is greatly appreciated. For their continued confidence and encouragement, the author is grateful to his family and friends. For her efficient and efficacious assistance with the technical aspects of producing this dissertation, the author is forever indebted to his wife, Margaret. Her patience, understanding and sustaining support have been salutary to making this work possible. iii
PAGE 4
TABLE OF CONTENTS Page ACKNOWLEDGMENTS i i LIST OF TABLES vii LIST OF FIGURES viii ABSTRACT xvi CHAPTER I A NEW PERSPECTIVE OF NONLINEAR ROUTING 1 Statement of the Problem 1 Statement of the New Perspective 3 Development of the New Perspective 5 Overview 5 Selection of a Primary Model 5 Calibration of the Primary Model to the Hydrologic System 6 Structure of the State Space Model 6 Calibration of the State Space Model to the Primary Model 7 The State Space Model as Surrogate for the Primary Model 8 An Application of the New Perspective 9 Organization of the Dissertation 11 II INSIGHTS INTO NONLINEAR ROUTING FROM LINEAR SYSTEMS THEORY 13 Linear Systems Theory 13 Introduction 13 Properties of a System..... 15 A General Linear Differential Equation 19 Transformation of Moments by a Stationary Linear System 22 Summary 2 5 Linearization of the SaintVenant Equations for a General Prismatic Channel 26 Insights from Linear Systems Theory 33 Introduction 33 The Inverse Gaussian PDF as an Analytical Solution of the Unsteady Flow Equations...... 33 Moments of the Inverse Gaussian PDF 34 The Inverse Gaussian PDF as an Inflow Hydrograph 37 Additional Properties of the Inverse Gaussian PDF 40 iv
PAGE 5
A Linear State Space Approach 41 Introduction 41 Criteria for Matching Impulse Response Functions 43 Computing the Parameters of a Gamma PDF 46 A Cascade of Linear Reservoirs State Space Model 52 Summary 57 III LIMITATIONS OF THE LINEAR THEORY 58 Introduction 58 Verification of the Linear SaintVenant Equations 59 Verfication of the Linear State Space Model for Small Inflow Transients 82 Comparison with the 3Parameter Gamma PDF 82 Comparison with the DWOPER Model 86 The Importance of Nonlinearities 98 Summary 1 1 IV A NONLINEAR STATE SPACE ROUTING MODEL 123 Analysis of the Important Nonlinearities 123 Nonlinear Characteristics of the Primary Model 125 Introduction 125 Theoretical Variations with Flow Level of the Impulse Response Function for the Linearized SaintVenant Equations for Incremental Inflow Transients 125 Simulated Variations with Flow Level of the Impulse Response Function for the Full SaintVenant Equations for Incremental Inflow Transients 130 Properties of the Nonlinear Storage Reservoir as a Building Block for the State Space Model... 137 Development of the Nonlinear Model 141 Introduction 141 The Final Configuration of the Nonlinear Model 148 The Path to the Final Model 153 Summary 158 V LIMITATIONS OF THE NONLINEAR STATE SPACE APPROACH 160 Introduction 160 Verification of the Nonlinear State Space Model.... 162 Introduction 162 SinglePeaked Inflow Hydrographs 162 A DoublePeaked Inflow Hydrograph 182 Continuity 182 The Number of Reservoirs 189 CPU Time 189 Summary 192
PAGE 6
VI HOW TO EMULATE ANY PRIMARY MODEL 193 Introduction 193 Select a Primary Model 193 Calibrate the Primary Model 193 Determine the Impulse Response Properties of the Primary Model 195 Estimate a, n and k as Functions of Flow Level 197 Parameters of the Nonlinear State Space Model 197 VII RELATED PREVIOUS WORK 200 Introduction 200 Related Routing Models 200 Observed Variations with Flow Level of the System Time Lag 203 Comparison with the MuskingumCunge Model 203 VIII CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 208 Conclusions 208 Recommendations for Future Work 211 APPENDICES A A REVIEW OF SOME HYDROLOGIC ROUTING MODELS 216 B A FORTRAN PROGRAM TO GENERATE HYDROGRAPHS WITH INVERSE GAUSSIAN PDF SHAPES IN A PRISMATIC CHANNEL 224 C A FORTRAN PROGRAM FOR THE NONLINEAR STATE SPACE ROUTING MODEL 227 BIBLIOGRAPHY 233 BIOGRAPHICAL SKETCH 239
PAGE 7
LIST OF TABLES Table Title Page, 3.1 Organization of Figures 3.1 through 3.18, Comparing the Linearized and Full SaintVenant Equations 61 3.2 CPU Times for the Numerical Solution of the SaintVenant Equations on a PRIME 750 Computer for Small Inflow Transients 83 3.3 Pairs of Figures with Identical Input Data Comparing the Linearized SaintVenant Equations and the Linear State Space Model with the Full SaintVenant Equations 90 3.4 Values of Parameters a, n and k for the State Space Model 91 3.5 CPU Times for the Linear State Space and DWOPER Models on a PRIME 750 Computer for Results of Figures 3.22 through 3.27 100 3.6 Organization of Figures 3.28 through 3.43, Comparing the Linear State Space Model with the Full SaintVenant Equations for Large Inflow Transients 102 5.1 Organization of Figures 5.2 through 5.17, Comparing the Nonlinear State Space Model with the Full SaintVenant Equations for Large Inflow Transients 163 5.2 Pairs of Figures with Identical Input Data Comparing the Linear and Nonlinear State Space Models with the Full SaintVenant Equations for Large Inflow Transients 164 5.3 Continuity Check for the Nonlinear State Space and DWOPER Models for Results of Figures 5.2 through 5.17 188 5.4 The Average Number of Reservoirs Used by the Nonlinear State Space Model in Figures 5.2 through 5.17 190 5.5 CPU Times for the Nonlinear State Space and DWOPER Models on a PRIME 750 Computer for Results of Figures 5.2 through 5.17 191 7.1 Values of m for Several Catchments 205 vii
PAGE 8
LIST OF FIGURES Figure 1.1 1.2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.1 3.2 3.3 3.4 Title Page An Overview of the New Perspective of Routing in Hydrology 4 An Overview of the Development of the New Perspective of Routing 10 A System ^ Calculation of Channel Outflow with the System Impulse Response Function 20 Development of the Analytical Impulse Response Function of the Linearized SaintVenant Equations 27 Channel Impulse Response Function Used to Represent Flow Hydrographs 3 5 Analysis of Limitations of Linear Theory 44 Comparison of the Inverse Gaussian and Gamma PDF with Parameters Matched by Method of Moments for a = 1, n = 4 and k = 2.0 48 Comparison of the Inverse Gaussian and Gamma PDF with Parameters Matched by Method of Moments for a = 8, n = 8 and k = 0.5 49 Comparison of the Inverse Gaussian and Gamma PDF with Parameters Matched by Method of Moments for a = 1, n = 1 and k = 0.5 50 Development of the Linear State Space Model 53 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q Q = 10000 cfs and a n = 3 hours 63 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for = 100000 cfs and a n = 3 hours 64 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 200000 cfs and a n = 3 hours 65 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 10000 cfs and a = 6 hours 66 in viii
PAGE 9
Figure Title Page 3.5 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 100000 cfs and a n = 6 hours 67 3.6 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 200000 cfs and o in = 6 hours 68 3.7 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 10000 cfs and a n = 12 hours 69 3.8 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 100000 cfs and a n = 12 hours 70 3.9 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for CL = 200000 cfs and a n = 12 hours 71 3.10 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for 0. = 10000 cfs and a n = 24 hours 72 3.11 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q Q = 100000 cfs and a n = 24 hours 73 3.12 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Rectangular Channel for Q = 200000 cfs and a ^ = 24 hours 74 3.13 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q = 10000 cfs and a n = 6 hours 76 3.14 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q Q = 100000 cfs and a in = 6 hours 77 3.15 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q Q = 200000 cfs and a in = 6 hours 78 3.16 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q = 10000 cfs and a n = 12 hours 79 3.17 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q = 100000 cfs and a n = 12 hours 80
PAGE 10
Figure Title Page 3.18 A Small Inflow Transient Routed with Linearized and Full SaintVenant Equations in a Triangular Channel for Q = 200000 cfs and a lQ = 12 hours 81 3.19 Verification Process for the Linear State Space Model 84 3.20 Comparison of the Linear State Space Model and 3Parameter Gamma PDF for a = 1, n = 4 and k = 2.0 87 3.21 Comparison of the Linear State Space Model and 3Parameter Gamma PDF for a 8, n = 8 and k = 0.5 88 3.22 A Small Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for in = 6 hours L = 10 miles and CL = 10000 cfs 92 3.23 A Small Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for in = 3 hours, L = 100 miles and CL = 100000 cfs 94 3.24 A Small Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations In a Rectangular Channel for in = 12 hours, L = 100 miles and 0^ = 100000 cfs 95 3.25 A Small Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a in = 24 hours > L = 10 miles and C^ = 200000 cfs 96 3.26 A Small Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a. = 6 hours, L = 100 miles and n = 100000 cfs 97 In 3.27 A Small Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Triangular Channel for in = 12 hours L = 10 miles and a = 200000 cfs 99 3.28 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for in = 3 hours and L = 10 miles 103
PAGE 11
Figure Title Page 3.29 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. =3 hours and L = 400 miles 104 in 3.30 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a, = 6 hours and L = 100 miles 105 in 3.31 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. =6 hours and L = 400 miles 107 in 3.32 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 12 hours and L = 100 miles 108 in 3.33 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 12 hours and L = 400 miles 109 in 3.34 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. =24 hours and L = 100 miles Ill in 3.35 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 24 hours and L = 400 miles 112 in 3.36 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a. =6 hours and L = 100 miles 113 in 3.37 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a. =6 hours and L = 400 miles 114 in 3.38 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a. = 12 hours and L = 100 miles 115 in .ii
PAGE 12
Figure Title Page 3.39 A Large Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Triangular Channel for u, = 12 hours and L = 400 miles 116 in 3.40 A DoublePeaked Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 100 miles 118 3.41 A DoublePeaked Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 400 miles and 0^ = 100000 cfs 119 3.42 A DoublePeaked Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 400 miles and 0^ = 10000 cfs 120 3.43 A DoublePeaked Inflow Transient Routed with the Linear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 400 miles and 0^ = 200000 cfs 121 4.1 Approach to the Nonlinear Analysis 124 4.2 Analysis of the Nonlinear Characteristics of the Primary Model 126 4.3 Variation with Flow of the Channel Lag for a. = 12 hours In for a Rectangular Channel 133 4.4 Variation with Flow of the Channel Lag for a. = 12 hours for a Triangular Channel 134 4.5 Variation with Flow of the Channel Standard Deviation for a = 12 hours for a Rectangular Channel 135 4.6 Variation with Flow of the Channel Standard Deviation for a n = 12 hours for a Triangular Channel 136 4.7 Variation with Flow of the Channel Coefficient of Variation for a n = 12 hours for a Rectangular Channel 138 jJL
PAGE 13
Figure Title Page 4.8 Variation with Flow of the Channel Coefficient of Variation for a = 12 hours for a Triangular Channel 139 4.9 Approach to the Development of the Nonlinear State Space Model 143 4.10 A Cascade of Nonlinear Reservoirs Governed by S m = ic0 144 4.11 Development of the Nonlinear State Space Model 147 4.12 A General Strategy to Emulate Any Primary Model 159 5.1 Verification Process for the Nonlinear State Space Model 161 5.2 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for o". =3 hours and L = 100 miles 165 in 5.3 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. =3 hours and L = 400 miles 168 in 5.4 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for in = 6 hours and L = 10 miles 17 5.5 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. =6 hours and L = 400 miles 17 1 in 5.6 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 12 hours and L = 100 miles 172 in 5.7 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 12 hours and L = 400 miles 173 in
PAGE 14
Figure Title 5.8 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 24 hours and L = 100 miles 175 in 5.9 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for a. = 24 hours and L = 400 miles 176 in 5.10 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a, = 6 hours and L = 100 miles 177 in 5.11 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a. = 6 hours and L = 400 miles 178 in 5.12 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a, = 12 hours and L = 100 miles 180 in 5.13 A Large Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Triangular Channel for a. =12 hours and L = 400 miles 181 in 5.14 A DoublePeaked Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 100 miles 183 5.15 A DoublePeaked Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 400 miles and 0^ = 100000 cfs 184 5.16 A DoublePeaked Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 400 miles and Q = 10000 cfs 185 5.17 A DoublePeaked Inflow Transient Routed with the Nonlinear State Space Model and the Full SaintVenant Equations in a Rectangular Channel for L = 400 miles and 0^ = 200000 cfs 186 6.1 How to Emulate Any Primary Model 194 vfji
PAGE 15
Figure Title Page 6.2 A Piecewise Constant Versus Flow Relationship for m 198 7.1 Lag Versus Mean Discharge Relationship, South Creek Experimental Catchment, University of New South Wales 204 7.2 A Comparison of the Nonlinear State Space and MuskingumCunge Models with the Full SaintVenant Equations at L = 400 miles in a Rectangular Channel... 207 vv
PAGE 16
Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A NEW PERSPECTIVE OF NONLINEAR ROUTING IN HYDROLOGY By George Francis Smith December 1983 Chairman: Wayne C. Huber Major Department: Environmental Engineering Sciences Existing routing models for hydrologic systems have either 1) a complex, nonlinear structure not amenable to second order (i.e., error) analysis, requiring considerable computational resources for solution; or 2) a simple, linear structure that is strictly valid over only a small flow range about which the model is linearized. The new perspective for routing surface runoff or channel flow uses any appropriate model, called the primary model, to simulate the observed physics for the catchment or channel of interest, and then emulates the mathematical behavior of the primary model with a computationally efficient nonlinear state space model that has structural properties amenable to second order analysis. The overall strategy matches impulse response properties of the state space model to those of the primary model over a wide range of flow levels. In general, impulse response properties vary with flow level. xv i
PAGE 17
Two important response properties are time lag and attenuation. Nonlinearities in time lag occur because water flows faster as it flows deeper. Nonlinearities in attenuation occur because deep waves disperse more rapidly than shallow waves. These two properties can be expressed as linear or piecewise linear functions of discharge in logspace, thus relating them directly to the parameters of a nonlinear state space model. The limitations of this approach to routing were examined by comparing results from the nonlinear state space model with a numerical solution of the SaintVenant equations. No significant differences were observed. The state space model required two to three orders of magnitude less CPU time than numerical solution of the SaintVenant equations A. generalized, computationally efficient state space routing model that can emulate the nonlinear behavior of any other routing model for downstream waves has been developed. Further studies are required to determine the full usefulness of this modelling approach. rvH
PAGE 18
CHAPTER I A NEW PERSPECTIVE OF NONLINEAR ROUTING Statement of the Problem Existing technology for routing unsteady flow in open channels and for routing surface runoff in a catchment has not proved adequate for some important hydrologic applications. Current routing methods are either inefficient when numerous simulations are required, as in flood frequency analysis; inadequate when the forecast requires an estimate of uncertainty in past or future flow conditions, as in water supply forecasting; or unsatisfactory when the information content of current data is important, as in realtime river forecasting. The above applications require 1. an objective estimate of uncertainty, 2. repetitive computations, and/or 3. consideration of second order (i.e., error) properties of models and data as well as the first order (i.e., mean) trajectory in time and space of the flow conditions; as well as requiring models or computational algorithms to account for significant nonlinearities of the channel or catchment response. For example, the National Weather Service (NWS), the Bureau of Reclamation, and the Federal Emergency Management Authority typically study the thousands of dams throughout the country by simulating various dam failure scenarios. Numerical solution of the full SaintVenant equations can impose a significant computational burden in such
PAGE 19
contingency dam break simulations for which repetitive calculations are required. Simpler models, while computationally efficient, do not properly account for the nonlinear behavior important to these analyses. The NWS studies of various possible dam failures are designed to provide timely and accurate forecasts of expected flow levels to communities downstream of the dam in the event of a dam failure. When the threat of a dam failure is imminent, the NWS Forecast Office can provide worst case (i.e., instantaneous failure) and alternative projected flow conditions to the media and emergency authorities. Several computer simulations involving flow routing downstream from each dam studied are required to produce that information. For daily, realtime river forecasting, as performed at NWS River Forecast Centers, timely processing of all current observations is critical. A rapid mechanism is needed to update the conditions projected by the computer model to account for the most recently observed data. A state space model of the full SaintVenant equations, which would be amenable to updating with modern estimation theory, imposes a large computational burden. Simpler routing models are not hydraulically adequate over the range of inflow transients encountered in flood forecasting. Although some simple models could be formulated for use with updating algorithms, their forecasting capability is diminished. by the limited flow range over which they properly model the nonlinear nature of water movement. Flood frequency studies, such as those performed during the hydraulic design of channel structures, often require computer simulations for many events over a wide range of flow conditions. Again, numerical solution of the full SaintVenant equations is costly
PAGE 20
and the simpler, linear models do not properly simulate for large variations of inflow. Statement of the New Perspective The new perspective of routing in hydrology is to let any model, called the primary model, simulate the physics of a channel or catchment, and to emulate the behavior of the primary model with a model structured to meet the applications described above. This process is depicted in Figure 1.1. To account for the modelling of second order properties (i.e., variance and covariance measurements of uncertainty) and to simplify computation, a state space form has been chosen as the model to emulate primary model behavior. The state space model structure is X =FX+GÂ„U, (I'D t+1 X,t t X,t t+1 where X^ = the system state vector, F v f = th e state transition matrix, A., t G^ t = the input coefficient matrix, and U r+1 = the vector of system inputs. The structure of a state space model is possibly the only one that meets the previously described objective of simplicity and has the potential for simulating second order properties. The state space model can be used with modern estimation theory in a Kalman filtering algorithm to propagate and update both mean and variance properties of the channel or catchment. The primary model simulates the physics of the hydrologic system while the state space model, because of its desired computational and structural properties, emulates the mathematical behavior of the primary model.
PAGE 21
REAL WORLD LINEAR SYSTEMS THEORY PRIMARY MODEL OF THE OBSERVED PHYSICS NONLINEAR ANALYSIS \ \ i / NONLINEAR STATE SPACE MODEL EMULATE THE MATHEMATICAL BEHAVIOR OF THE PRIMARY MODEL Figure 1.1 An Overview of the New Perspective of Nonlinear Routing in Hydrology
PAGE 22
The essence of this new perspective of routing is its twostep approach to the problems encountered with existing routing technology. First, the physics of the system are modelled with a primary model, then the mathematics of the primary model are emulated by a state space model. The state space model is a model of the primary model. Development of the New Perspective Overview Several steps are required to apply this new perspective to a routing problem of interest. The steps are 1. select a primary model, 2. calibrate the primary model to the hydrologic system, 3. determine the appropriate structure for the state space model, 4. calibrate the state space model to the primary model, and 5. emulate the primary model behavior with the state space model. Selection of a Primary Model The primary model simulates the physics of a channel or catchment. Because, in general, the forces governing the motion of water interact in a complicated manner, a complex model typically will be most appropriate to describe the physical properties of the hydrologic system. However, the primary model does not have to be complex. Any simple model can be chosen to simulate the system physics based on data or computational limitations, or modeller preference. A
PAGE 23
summary of some available channel routing models is presented in Appendix A. Calibration of the Primary Model to the Hydrologic System The effort required to calibrate the primary model will vary greatly depending upon the complexity of the primary model chosen. If a physically detailed hydraulic model is selected, channel crosssection data must be obtained and reduced to the form required by the model. This requires a complicated process of adjusting friction factors, determining active versus inactive flow area, etc. For simpler models, the calibration process usually requires determining a few parameters based on the comparison of simulated and observed data for some past inflows. Parameter selection for simple models can often be automated, given inflow and outflow records. Even for some complex hydraulic models, final parameter adjustment can be performed automatically (Fread and Smith 1978). Structure of the State Space Model The mathematical properties of the primary model can be represented in numerous ways, given the general structure of equation 1.1. The system states (i.e., the elements of the vector X_) can represent many possible quantities. The states of the surrogate model are determined by the important relationships in the primary model. State space models of hydrologic systems can be constructed with ^ physically based states such as the quantity of water in a reservoir (SzollosiNagy 1981) or with states which are purely mathematical
PAGE 24
artifacts that reproduce observed system behavior (Goldstein and Larimore 1980). In fact, the full SaintVenant equations can be expressed in state space form. At each time step, values of flow and crosssectional area at numerous points along the channel represent the state of the system. The state vector, Xj. could be constructed from these flow and area values. However, since the number of points defined along a channel typically is large, the state vector would be correspondingly large. For example, some NWS applications of the Dynamic Wave Operational Program (DWOPER), which numerically solves the SaintVenant equations (Fread and Smith 1978) have 25 and 45 points defined along the channels, respectively, for the ColumbiaWillamette River system and the OhioMississippi River junction. The state vectors for these systems would consist of 50 and 90 elements, respectively, with state transition matrices of 2500 and 8100 elements. Obvious dimensionality problems arise in the computer solution of equation 1.1 for systems of these sizes. The structure of the nonlinear state space model developed is presented in detail in Chapter IV. The approach taken uses the nonstationary linear structure of the state space model as a surrogate for the stationary nonlinear behavior of the primary model. Calibration of the State Space Model to the Primary Model Problems encountered in certain routing applications can be solved by using the state space model to emulate the mathematical behavior of the primary model. Since the state space model is calibrated to the mathematical properties of the primary model, many approaches are possible.
PAGE 25
The approach taken here draws on the general, powerful theory of linear systems. As described in Chapter II, the nearly linear behavior of the primary model for a small flow range is emulated with an equivalent linear system. This is accomplished by matching the impulse responses of the primary and state space models for small flow ranges about various reference discharges. The inherent nonlinearities that cause the behavior of the state space model to agree with or to depart from the primary model behavior are variations with discharge of the lag and dispersive properties of the impulse response function. These properties vary substantially over large changes of inflow. A simple functional relationship can be derived that specifies the variation of F_ and ^with discharge to emulate properly the lag and dispersive properties. Given such a functional relationship, the state space model can emulate the full nonlinear and dynamic behavior of the primary model over all flow ranges. The State Space Model as Surrogate for the Primary Model Simple linear models have been used in place of complex routing models for many years. They replace the complex model, but do not properly represent some of the important properties of the complex model. With the new perspective, the primary model is not replaced. Instead its mathematical structure is emulated by a state space model with computational and structural properties appropriate for the solution of problems encountered in many routing applications. The state space model acts as a surrogate for the primary model.
PAGE 26
The state space model typically has much lower dimensionality than the primary model and imposes a much smaller computational burden. Because the state space model emulates the full nonlinear behavior of the primary model, it can be used in applications where the primary model is not appropriate for error analysis because of the computational burden or improper structure. This approach presents a model mathematically equivalent to the primary model, that can be used on smaller computers than would be possible with the primary model itself. An Application of the New Perspective To explore the new perspective's applicability to routing in hydrology, the general steps described above were followed for a specific primary model. As shown in Figure 1.2, the application selected was the full, nonlinear SaintVenant equations as the exemplar primary model to simulate downstream flow in a prismatic channel. This primary model and physical system were chosen because 1. the SaintVenant equations are generally recognized as the most complete representation of unsteady flow phenomena, 2. an analytical solution for the impulse response of the SaintVenant equations can be derived for downstream wave movement in prismatic channels, and 3. numerical solution of the SaintVenant equations is wellstudied, and many excellent solution algorithms are available for comparison with the surrogate model (Liggett and Cunge 1975). Analysis of the ability of a state space model to emulate the behavior of the full, nonlinear SaintVenant equations should lend
PAGE 27
10 NEH PERSPECTIVE QF ROUTING I SELECT SFIINTVENANT EQUATIONS AS EXEMPLAR PRIMARY MODEL LINEAR ANALYSIS OF SAINTVENANT EQUATIONS NONLINEAR ANALYSIS OF SAINTVENANT EQUATIONS DEVELOP LINEAR STATE SPACE MODEL i VERIFICATION QF LINEAR STATE SPACE MODEL DEVELOP NONLINEAR STATE SPACE MODEL VERIFICATION QF NONLINEAR STATE SPACE MODEL \ y GENERAL PROCEDURE TO EMULATE ANY PRIMARY MODEL Figure 1.2 An Overview of the Development of the New Perspective of Routing
PAGE 28
11 insight into the applicability of this approach for other primary models. Because, in general, no analytical solution for the impulse response of the primary model is available, the case chosen allows comparison of the true impulse response with results obtained numerically. Organization of the Dissertation The organization of this dissertation is schematically presented in Figure 1.2. The new perspective of routing and the reasons for selecting the SaintVenant equations as the exemplar primary model have been presented. A linear analysis of the SaintVenant equations and development of a linear state space model follow in Chapter II. In Chapter III, results from the linear state space model are compared with the full SaintVenant equations for small inflow transients about a reference flow level. These results indicate that the linear state space model is a very adequate surrogate for the SaintVenant equations over small flow ranges. Further nonlinear studies are motivated by comparison of the linear state space model with the primary model for large inflow hydrographs. Nonlinear behavior of the primary model causes the comparison of results with the linear state space model to deteriorate. A nonlinear analysis of the SaintVenant equations primary model and the development of a nonlinear state space model are conducted in Chapter IV. Chapter V is devoted to the verification of results for the nonlinear state space model versus the primary model for large inflow hydrographs. Results from the primary and surrogate models for a complex inflow hydrograph were comparable. However, the nonlinear state
PAGE 29
12 space model used two to three orders of magnitude less CPU time than the numerical solution of the SaintVenant equations. A general procedure to emulate any primary model with the nonlinear state space model is presented in Chapter VI. By following the steps here described, the nonlinear state space model can be used as a surrogate for any channel routing or surface runoff model. The relationship of the current work to previous studies of routing in hydrology is the subject of Chapter VII. A review of pertinent work, by other authors, some data to support the power function relationship between flow level and time lag, and a comparison of the nonlinear state space and MuskingumCunge models with the full SaintVenant equations are the major topics of this chapter. Conclusions on the new perspective of routing in general, and on its application to the exemplar primary model in particular, are stated in Chapter VIII. Recommendations for future work, Including the necessary tasks for the application of estimation theory with the nonlinear state space model and possible alternative approaches for representing the nonlinear model structure, form the final portion of this last chapter.
PAGE 30
CHAPTER II INSIGHTS INTO NONLINEAR ROUTING FROM LINEAR SYSTEMS THEORY Linear Systems Theory Introduction Linear systems theory provides many tools for the analysis of linear phenomona. Linear systems theory, strictly applied, would have few uses because, in general, the only truly linear systems are mathematical abstractions. However, the spirit of linear systems analysis lends itself to unraveling those aspects of nonlinear systems that behave in a nearly linear manner. The phenomenon of routing in hydrology is, in general, a very nonlinear process. Nevertheless, important insights into hydrologic system behavior may be gained by applying some of the tools of linear systems theory. The concepts of linear systems theory have been applied to many and varied systems and are well known. Those concepts germane to this work are summarized below. Systems can be thought of abstractly as the mechanisms that interrelate two objects (Dooge 1973, pp. 34). The operational definition used in this work is that a system is the process that transforms input signals into output signals. For the hydrologic processes analyzed here, the complex relationships among the forces governing the motion of particles of water are conceptualized in a system as shown in Figure 2.1. The value of this view of a system lies 13
PAGE 31
14 INPUT SIGNALS
PAGE 32
15 in the analytical tools it provides since the mathematical properties of the system determine which tools are applicable. Properties of a System The properties of consequence for hydrologic systems are whether the system is a. determinsitic or probabilistic, b. distributed input/distributed output c. causal, d. stationary, and e. linear. Typically the response of a timevarying, nonstationary, distributed, probabilistic, nonlinear hydrologic system is simulated by a time invariant, stationary, lumped, deterministic, linear model. And we wonder why there are discrepancies between observations and simulated outputs! A system is determinsitic if the precise state of the system (the values of X in equation 1.1) can be foretold. If a random component is introduced into the variation of the system states, the system becomes probabilistic. The inputs and outputs of most hydrologic systems do not interact at a single point. As such, there should properly be some spatial description of the system processes. The complex, dynamic, nonlinear models of channel routing of surface runoff typically account for some of the spatial variation of parameters over the channel or catchment. The parameters of simple linear models usually apply over the entire channel or catchment area and, therefore, such models operate as if all
PAGE 33
16 input/output interactions occurred at a single point (i.e., they are lumped models) In general, a system can have several inputs and outputs as shown in Figure 2.1. The hydrologic systems analyzed in this work are assumed to have a single input and a single output. A causal system is one in which outputs cannot precede inputs. All physical systems are causal. This restriction limits the range over which integral equations describing the system are applied. The stationarity of a system is defined by whether or not the input/output relationship changes with time. If a system is stationary, identical input produces identical output. Linearity is the property that opens the doors to the storehouse of mathematical tools. If system inputs are defined as x lf system outputs as y t > and the transformation of x into y by the system is symbolized as > then given _> y and x > y 11 2 2 a system is said to be linear if a x + a Â• x > a Â• y + a Â• y ( 2>1 ) 11 2 2 112 2 where a and a 2 are constants (Thomas 1969, p. 140). The analysis of nonlinear systems is much more difficult than that of linear systems because the principle of superposition applies to linear systems (Thomas 1969, p. 139). Superposition, which follows directly from the definition of linearity, says that if inputs are added or scaled, outputs are similarly added or scaled. Therefore, the system behavior can be analyzed independently of the inputs. This is a point of great significance: if the system output for an instantaneous unit
PAGE 34
17 input can be found, the system behavior for any_ input can be obtained by considering the actual input as a series of scaled unit inputs at various instants of time (Liu and Liu 1975, pp. 135146). Given the system output for an instantaneous unit input (called the impulse response), the output for a stationary linear system can be obtained for any input with the convolution integral ,(t) = / h(tT)x(T) dT (2.2) where x('), y(* ) are the system inputs and outputs respectively, and h( ) is the system impulse response (Johnson and Johnson 1975, pp. 137140). Stationary nonlinear behavior can sometimes be simulated by a nonstationary linear model. This is, in fact, exactly the strategy behind the nonlinear state space model developed in Chapter IV. The theoretical basis is most clearly seen by recognizing that for a linear model, the input/output relationship can be represented by matrix multiplication as I ] 1,1,1, 12 3 h h h 11 12 13 h h 22 23 h 33 1" (2.3) = [o , Â• Â• Â• o J 12 3 n where I i = the input at time i, an element of row vector U = the output at time j, an element of row vector 0, and h j = an element of the input/output relationship matrix _H.
PAGE 35
18 Note that I and are time vectors, in contrast to the state vector X. of equation 1.1. Matrix _H transforms I_ into 0. What we want at any time, ..Â„ n in vector 0. To do this we need the t, is to generate one entry, t LU veuut 1L' appropriate weights to apply to all inputs 1^ U p to and including l t These weights will fill one column of matrix H. However, it is physically impossible to observe the weighting function of a system. An approach which circumvents that problem is to generate the model impulse response. This is equivalent to solving f or 0_ in equation 2.3 when I_ is all zeroes except for a single entry of 1 in, say, the 1 th location. Call this vector 6_ The output, 0_, generated will be the i c _row of matrix H. The input, 6^ selects the ltb row of H and presents it as the output, 0. For a linear nonstationary system the rows of H_ will, in general, be different from each other (i.e., each row, i, will represent the system impulse response at time i). The assumption implicit when using a single system impulse response in the convolution integral is that the system is linear and stationary. In this case, the i th row of H will be equal to the (il) th row of H_ shifted one position to the right. Therefore, equation 2.3 can be rewritten for a stationary linear system as h 1
PAGE 36
19 The desired weighting functions (i.e., the columns of H) are seen in equation 2.4 to be the timereversed images of the system impulse response for a stationary linear system. Figure 2.2 shows how the convolution integral (equation 2.2) exploits this to compute the system output. Nonlinear system behavior can be modelled by letting the columns of H in equation 2.3 vary each time step. The elements on any diagonal will not be equal, as is the case for the linear stationary system represented by equation 2.4. Therefore, the system impulse response (the rows of H) will change each time step, thus permitting a nonstationary linear model to account for stationary nonlinear behavior. The state space model developed in Chapter IV emulates the nonlinear system behavior by varying the weighting function based on the level of flow. A General Linear Differential Equation Linear timeinvariant equations can be used to approximate the behavior of many physical systems. This is a good approximation when system characteristics change very slowly relative to variations of the inputs (Ogata 1967, p. 307). Channel and catchment systems can be well represented by linear timeinvariant equations. A general differential equation for a nonstationary linear system can be written as ( t ).^l2 + a l( t).il^l + Â• Â• Â• n n1 j,111 dt dt (2.5) + ai (t).l4r+ a Q (t)y(t) = x(t)
PAGE 37
20 fri o TIMEREVERSED IMAGE OF THE SYSTEM IMPULSE RESPONSE fcÂ— (t)0 TIME INTO THE PAST (} ttc)t THE CHANNEL OUTFLOW AT ANY TIME IS A WEIGHTED AVERAGE OF THE PAST AND ; PRESENT INFLOWS. THE TIMEREVERSED I IMAGE OF THE SYSTEM IMPULSE RESPONSE ; REPRESENTS THE WEIGHTING FUNCTION. v> 2 TS, (PRESENT] N 0(t)< /hC<)'I(t<)ck V (PRESENT) Figure 2.2 Calculation of Channel Outflow with the System Impulse Response Function
PAGE 38
21 where x(t) = the system input, y(t) = the system output, a (t), ai (t), . ., ani(t), a^O are coefficients, and the initial conditions are specified for .n1. y( O, &&K .... and Â•y(t ) dt dt n1 The system is linear because there are no powers or products of y. Equation 2.5 can be rewritten in matrix form as djY(t), dt A(t)Y(t) = B(t)X(t) (2.6) where Y(t) = y(t) dy(t) dt d n "Vt) dt n1 A(t) = a (t) "aTO a
PAGE 39
22 X(t) = x(t) If the system is stationary the coefficients A_(t) and _B(t) of the underlying differential equation 2.5 are not functions of time. In this case equation 2.6 can be rewritten as fLlKt)) A Y(t) = B x(t) (2.7) dt which can be solved using the integrating factor e to yield Y( t ) =e (t t0) ^Y(t ) + J e (t T) ^B.X(x) dx (2.8) (Ogata 1967, p. 315). The first term on the right hand side of equation 2.8 is the zero input response of the system. With no input the system would move toward equilibrium starting from this point when t = t Q The second term is the zero state response, which describes the motion of the system driven by input X(x ) for t Q <^t <_ t Transformation of Moments by a Stationary Linear System One way to describe a system's input and output signals is to study their moments. The following section, on the use of transform methods to study moments of stationary linear systems, is based on work by Dooge (1973, pp. 113115). The R th moment of a function about the time origin is
PAGE 40
23 u:(f) = / f(t)t R dt (2.9) and the R th moment about the center of mass Is U n (f) = / f(t)(tU') R dt (2.10) R oo 1 The output signal of a stationary linear system can be expressed with the convolution integral (equation 2.2) in terms of the input signal and the system impulse response. The relationships between moments of the input, output and impulse response are most easily seen by transforming equation 2.2 with either the Fourier or Laplace transform. Since moments about the center of mass of the input and output signals are of interest, we will use the bilateral Laplace transform defined by 00 F D (s) = / f(t)e" st dt (2.1D _oo where f(t) is a function in the time domain. Differentiating equation 2.11 R times and setting s = produces the following expression for the R th moment about the origin in terms of the Laplace transform (Dooge 1973, p. 114). Â„.<Â„ h, 1 .A,w!h (7 l2) ds The convolution integral can be expressed in the transform domain by taking the bilateral Laplace transform of equation 2.2 to yield Y B (s) = H B (s)X B (s) ( 2 13 > where Y B (s), H B (s) and X B (s) are the Laplace transforms of y(t), h(t) and x(t), respectively (Dooge 1973, p. 114).
PAGE 41
24 The R th moment of the output about the origin follows from equation 2.12 Â„. (y) c R .4[V>U (2 14) ds Substituting equation 2.13 for Y B ( S ) gives U'Cy) = (i) R 4t H B (8) X B (8) ls0 (2 15) ds Solving equation 2.15 with Leibnitz's rule for continued differentiation of a product, and substituting for the moments with equation 2.12 yield R k=0 where f R ) indicates a combination of R things taken k at a time, or k y r R l = R! l k J k!(Rk)! U' (h) = the (Rk) th moment about the origin of the system Rk impulse response, and U'(x) = the k th moment about the origin of the system input k signal (Dooge 1973, p. 114). If, instead of about the origin, the moments are taken about the centers of mass of the input, output and system impulse response functions, the following relationship holds (Dooge 1973, p. 115) k=0 Only the first moment about the origin and the second moment about the center of mass are needed for the ensuing work. For R = 1 equation 2.16 reduces to U'(y) = U'(h) + U'(x) (2.18) 1 1 1
PAGE 42
25 Since U (Â• ) = 0, for R = 2 equation 2.17 becomes U (y) = U (h) + U (x) (2.19) 2 2 2 Therefore, for the first moment about the origin (the time lag) and the second moment about the center of mass (the variance), the moment of the system output is the sum of the moments of the system impulse response and the system input (Dooge 1973, p. 115). Summary The impulse response is the output generated by an instantaneous unit input to a system. Because the impulse response function completely describes the behavior of a stationary linear system, we would like to determine the impulse response of the channel or catchment physical systems. The best description of the real world that we have is the primary model. The goal of linear analysis of the system is to determine the impulse response of the primary model at various reference flow levels. The changing properties of the impulse response function with flow level lend insight into the nonlinear behavior of the physical system. An operational way to determine the impulse response of the primary model is to pulse the model with a unit input and observe the output. That is, in general, the approach that will be taken to find a system's impulse response. An analytical solution can be found for the impulse response of the SaintVenant equations of unsteady flow for downstream waves in a prismatic channel. To derive the impulse response, the SaintVenant equations must be linearized about a reference discharge. This process is described in the next section.
PAGE 43
26 Linearization of the SalntVenant Equations for a General Prismatic Channel To determine analytically the impulse response for routing unsteady flow, the SaintVenant equations must be linearized in terms of small perturbations about a reference flow condition. The steps presented in this section are outlined in Figure 2.3. Dooge (1973, pp. 245253) presented the special case for linearizing the SaintVenant equations for a rectangular channel and Chezy friction formula. The ensuing derivation follows Schaake (1980) who extended Dooge's work for a general prismatic channel and resistance formula. The SaintVenant equations for unsteady flow can be expressed as d A + ll= Bll+i^lq (220) 3t
PAGE 44
27 LINEARIZE ABOUT FLOW LEVEL, % REPRESENT CHANNEL GEOMETRY AS B 8 *A r REPLACE HIGH ORDER TIME DERIVATIVES WITH SPACE DERIVATIVES USING KINEMATIC WAVE CELERITY SAINTVENANT EQUATIONS i LINEAR, SECOND ORDER HYPERBOLIC EOUATON LINEAR, SECOND ORDER PARABOLIC EQUATION IMPULSE RESPONSE IS INVERSE GAUSSIAN PDF Figure 2.3 Development of the Analytical Impulse Response Function of the Linearized SaintVenant Equations
PAGE 45
28 <; = the channel bottom slope (Henderson 1966, p. 287). Equation 2.20 is the continuity equation for a channel. The effects of the forces acting on a control volume within the channel are modelled by equation 2.21, the momentum equation. Since the continuity equation is already linear, only the momentum equation requires linearization. Equation 2.21 is first rewritten in terms of the dependent variables Q and A. To this end, the friction term and channel geometry are expressed as v = Q=k.() P ./s7 (2.22) and B = BA r (2.23) where B is the top width, B Q is the top width for flow Qq and k, p and r are described below. Equation 2.22 is a general relationship for friction in open channel flow which can be related to the wellknown Manning or Chezy friction formula. With the approximation R A Â• A (2.24) R h P B where R^ = the hydraulic radius and P = the wetted perimeter. Equation 2.22 becomes Chezy' s friction formula with k = C a (2.25) and p = 1/2 or Manning's friction formula with
PAGE 46
29 k = 1.49/n and < 2 26 > p = 2/3 where C = the Chezy flow resistance factor and n = the Manning flow resistence factor (Daily and Harleman 1966, pp. 297298). Equation 2.23 can approximate most regular cross sectional shapes. For a rectangular channel r = 0; while for a triangular channel r = 1/2. However, a trapezoidal channel is not well modelled by equation 2.23 for all flow ranges. For shallow flow a trapezoidal channel behaves similarly to a rectangular channel; while for deep flow it can be approximated by a triangular channel. There is a transition flow range where equation 2.23 does not represent the change in top width with area for a trapezoidal channel. Substituting 2.23 into 2.22 and rearranging gives = 5 f k 2 A 2lp(lr)+lJ (2.27) Let m = p(lr)+l and Ws~ (2.28) (2.29) where m is analogous to the dimensionless kinematic wave parameter defined by Eagleson (1970, p. 250). Substitution of 2.28 and 2.29 into 2.27 gives
PAGE 47
30 (2.30) Substituting equations 2.20 and 2.30 into the momentum equation gives (g A W,.^ +2 B Q A .Â£ tB A ^ = (2.31) 3+r g'B 2 3+r2m g bA Â• S => Q Â• A ct 2 where henceforth the terms involving q will be ignored because we are studying prismatic channels with no lateral inflows. Finally, differentiating 2.20 with respect to x, differentiating 2.31 with respect to t while holding the lefthand coefficients constant, and equating the common ,, \% term, results in A jn 3 2 Q Â„ 3 2 Q 3 2 Q Â— v* J Â• 2* v 5 Â— 5Â— = 1 ax 2 3x3t 3t 2 (2.32) Â„ 3Q 2 g S 3Q 2mg'S ^+ ^^t 6 J i x v 3 1 Linearizing equation 2.32 about Q A, v Q and B Q forms the linear, second order, hyperbolic equation gA o^ 3 2 Q 3 2 3 2 Q B j 9x 2 3x3t 3t 2 (2.33) ?Â• ?Â•
PAGE 48
31 The impulse response of the linearized form of equation 2.33, presented in Appendix A, is a complex function containing a Dirac delta and a modified Bessel function (Harley 1967, p. 18). A simpler, approximate solution can be found with an assumption from kinematic wave theory. The 3x3t and 3 t 2 terms of equation 2.33 can be replaced by a Sx 2 term using the relationship 1Q _ r 30 (2.34) 3t 3~x where dx c = 1 Â— (the wave celerity), dt This substitution is routinely made in the derivation of diffusion routing equations (Koussis 1976). Here its use results in the simplified linear approximation to the complete equations 9Q .3Q Q h^.BJO (2.35) 3t C 5 2.B q .S o 9x2 where $ 2 = F 2 [lm(2m)l (2,36) in which Q /A F = (the Froude number at Q ) (2.37) /gA /B and c = mv (the kinematic wave velocity) (Schaake 1980). (2.38) The diffusion form of equation 2.35 appears often in channel and pollutant transport modelling (Fischer 1967, Carslaw and Jaeger 1959, Bansal 1971). Koussis (1976) followed a derivation similar to the one
PAGE 49
32 presented above, but was interested in a simpler numerical solution to the SaintVenant equations, not an analytical one. The interest in this section is in an analytical solution to equation 2.35. This will provide the theoretical basis for comparison of results from a numerical solution of the full SaintVenant equations, and the linear and nonlinear state space models developed later in this chapter and in Chapter IV, respectively. The influence of all terms in the SaintVenant equations has been retained in the parabolic form of equation 2.35. The impulse response function of equation 2.35 is (xc t) 2 h(x t) = X > ex P {TT&} (2.39) 4ifD t 3 where D^_flS 2 ) (2' 4 ) 2'B S l ; (Harley 1966, Schaake 1980). The impulse response, h(x,t), defines the flow for all positive x and t. The area under h(x,t) is unity, (i.e., / h(x,t) dt = 1) to preserve continuity, and the units are time 1 Equation 2.39, the solution for an instantaneous input at x = 0, t = in equation 2.33, appears in many scientific disciplines as a solution to the advectiondif fusion equation (Carslaw and Jaeger 1959, p. 50). Since the diffusivity, 0, is constant for any reference flow, QÂ„ Fickian diffusion is represented by equation 2.39 (Slade 1968, pp. 8081).
PAGE 50
33 Insights from Linear Systems Theory Introduction The analytically derived impulse response for the SaintVenant equations is strictly correct only for small variations about the reference flow. In the spirit of linear approximations for nonlinear systems, the impulse response can be used to approximate the channel response for all flow ranges, as is shown in Chapter III. Equation 2.39 completely describes the response of the general prismatic channel for flows close to Q This analytic solution can be manipulated to gain insights into the behavior of the SaintVenant equations, from which 2.39 was derived. The Gaussian form of the solution of equation 2.39 for fixed t is well known (Crank 1956, pp. 911). However, Schaake (1980) recognized that for fixed x the form of equation 2.39 is an inverse Gaussian probability density function (pdf ) Given this insight, a wealth of information can be obtained from previous studies of this statistical distribution (Johnson and Kotz 1970, pp. 137153). The Inverse Gaussian PDF as an Analytical Solution o f the Unsteady Flow Equations Since 1871, when Barre De SaintVenant proposed the equations of unsteady flow in an open channel, people have sought an analytical solution for equations 2.20 and 2.21. No analytical solution to the full SaintVenant equations has been found. An analytical expression for the outflow hydrograph would be ideal for studying the behavior of a channel. However, the ideal has not been attained because, in addition
PAGE 51
34 to the lack of an analytical solution of the SaintVenant equations, there is no analytical expression for a general inflow hydrograph. There is, however, an analytical solution to the equations of unsteady flow for a special class of inflow hydrographs in a prismatic channel. Because equation 2.39, the equation for an inverse Gaussian pdf, is the response for an impulse input to a prismatic channel, an inverse Gaussian pdf can be observed at any point along the channel. Parameters of the inverse Gaussian pdf are determined by the channel properties and the distance, x, from the point at which the impulse is input. As shown in Figure 2.4, an inverse Gaussian pdf with x = x Q + l is the analytical expression for the outflow hydrograph when an inflow hydrograph is specified by an inverse Gaussian pdf with x = x Q Channel properties determine the parameters Cq and D of equation 2.39. This means that for the class of inflow hydrographs defined by inverse Gaussian pdfs, the outflow hydrographs from a prismatic channel are also inverse Gaussian pdfs (Schaake 1980)! This is of great significance because the parameters of the outflow inverse Gaussian pdf are related simply to the parameters of the inflow inverse Gaussian pdf and the channel properties. Moments of the Inverse Gaussian PDF The time delay and dispersion characteristics of 2.39 are expressed by its first two moments (Johnson and Kotz 1970, pp. 137140). The first moment about the origin (i.e., the mean) describes the time lag properties as k = t. = i I c (2.41)
PAGE 52
35 A\ caW in in 5 gfe x /K CO r coÂ£
PAGE 53
36 The second moment about the mean (i.e., the variance) describes the dispersion properties as Q x k = a 2 = S >(l*' 2 B Â• S Â• C 3 (2.42) The coefficient of variation, a dimensionless measure of the impulse response, can be obtained from 2.41 and 2.42 as 1/2 B Â• S Â• c x 14' (2.43) The form of the coefficient of variation, equation 2.43, indicates that the kinematic wave assumption is most nearly correct for shallow U flow (small < ) on steep slopes (large S Q ) Woolhiser and Liggett Q ,. ..uunu,,,, .x,,.Â„ ^gget (1967) define a dimensionless parameter to determine when the kinematic wave assumption is appropriate (2.44) L S k = _0 0_ WL H .p2 where L = characteristic length, S n = channel bottom slope, Hq = depth of flow, and F 2 = Froude number. If k WL > 10, the SaintVenant equations are approximated well by a kinematic wave solution (Eagleson 1970, p. 336). These conditions are generally met by overland flow, with its shallow depth of flow and steep slopes. For river channels with steep slopes, the S Q term of the momentum equation dominates, making the kinematic wave assumption a good
PAGE 54
37 approximation of the SaintVenant equations. For rivers with flat slopes, the S and ^terms are dominant in equation 2.21 (Henderson 1966, p. 364). The diffusion analogy model, discussed in Appendix A, uses only these two terms for channel routing. The Inverse Gaussian PDF as an Inflow Hydrograph The interrelationships between peak flow (Q_), celerity (c Q ) and timing properties (as measured by the variance, o 2 ^) of the inflow hydrograph can be expressed for a given total flow volume (V) (Schaake 1980). With these relationships several properties of hydrographs input to routing models can be controlled. Any three of Q p i c a 2 or V can be specified arbitrarilv. Since equation 2.39 is the in solution for a pulse input at x = and t = 0, we will define a point downstream at x = Xq to observe the flow for t >_ 0. With this notation and h(,) defined by equation 2.39, the total flow volume, Inflow hydrograph, and peak discharge are, respectively, V = f Q(x ,t) dt = / Vh(x ,t) dt = VJ h(x ,t) dt o (2.45) Q(x ,t) = Vh(x ,t) ( 2 46 ) o Q = Vhfx ,t (x ), y p v P where h(x,t) is a maximum when (2.47) t p (,>iÂ£i{i*^C) 2 i'" > i 3D ; c Â—
PAGE 55
38 By choosing any three of V, Q p a\ n and Cq a hydrograph with certain desired properties can be selected. Since equation 2.39 defines the response to an impulse input, an inverse Gaussian shape hydrograph can be observed at any point along the channel. The channel can be divided into two sections as shown in Figure 2.4. The upstream, ficticious section has an impulse input and produces an inverse Gaussian pdf outflow hydrograph which becomes the inflow hydrograph for the downstream, real channel of interest. When an impulse is input to the upstream section at x = 0, an inverse Gaussian pdf appears as input to the downstream section at x = Xq The variance of the flow transient as x = x Q (a 2 ) can be selected by proper choice of Xq Letting a 2 and x of equation 2.42 equal a 2 and x Q respectively, gives a 2 Q X h^ 2 ] (2.49) Substituting equation 2.40 into equation 2.49 simplifies the expression for the inflow variance to (2.50) The parameters D and c. are specified by the channel geometry and slope. The value of xq can be arbitrarily selected to give the desired inflow variance. Equation 2.50 can be rearranged to solve for x Q as c 3 o 2 = _0 In (2.51) 2D 2
PAGE 56
39 The channel section where <_ x <_ Kq is a ficticious reach designed solely to generate an inflow hydrograph with desired properties for the real channel being modelled. This real channel is defined by the x _<_ x <_x Q + L channel section where L is the length of the real channel. The outflow hydrograph is also an inverse Gaussian shape with variance 2'D'fx + L] a' _0 (2.52) out c 3 Equations 2.18 and 2.19 can be rewritten to solve for the mean and variance of the impulse response U'(h) = U'(y) U'(x) (2.53) 1 1 1 and U (h) = U (y) U (x) (2.54) 2 2 2 With these relationships and equations 2.41 and 2.50, the mean of the impulse response function for the linearized SaintVenant equations can be expressed as x + L x u(h) =_Q JL (2.55) 1 c c 1 U'(h) = lag K i(2.56) i channel c 1 an d the variance of the impulse response function is 2D'(x + L) 2D'x 0(h) o2 o] B. 2. (2.57) 2 out in c 3 c 3
PAGE 57
40 Â„ ( h ) *Â£*. (2.58) o channel 3 Equations 2.56 and 2.58 express the lag and variance, respectively, which are added to an input hydrograph as it passes through a prismatic channel specified by parameters D, c Q and L. The channel impulse response will vary depending on the reference flow level chosen, since D and c are functions of Q Note, however, that the channel impulse response is independent of x Q This means that once Qq is chosen, the impulse response of the system is not dependent on the input signal. The development of a linear state space model which follows later in this chapter uses the properties of the linearized SaintVenant equations described above. The remainder of this section follows Schaake (1980) and focuses on several properties of equation 2.39 that lend additional insight to routing in hydrology. Additional Properties of the Inverse Gaussian PDF Attenuation of the peak inflow, Q to the real section of the channel can also be determined from an analysis of equation 2.39. Define the attenuation of the peak inflow at a point distance L below the upstream location, Xq as V'hfx + L,t (x + L)' A(L) = V'hfx ,t (x ), P g (2.59) where the numerator is the peak outflow rate at location x Q + L and the denominator is p Cancelling the common factor produces hfx + L,t (x + L)) A(L) = 3 L_J (2.60) h(x t (x )) P
PAGE 58
41 Schaake (1980) notes the interesting result that although A, the inflow hydrograph attenuation, is independent of V, the downstream peak outflow rate is directly proportional to V. The channel lag of the center of mass (i.e., the mean) of the inflow hydrograph is given by equation 2.56. The peak lag time can be derived from equation 2.39 as C Â„ = t fx + LJ t (x P* P^ P 1/2 1/2 'pi 1 + c (x + L 3^D + J "snr (2.61) (2.62) A Linear State Space Approach Introduction The inverse Gaussian pdf is a linear approximation of the full nonlinear SatntVenant equations. For small variations in flow the Inverse Gaussian pdf produces results nearly Identical with the full equations (Chapter III). Unfortunately, derivation of a state space model for equation 2.39 is cumbersome. A state space model with parameters related to the parameters of the inverse Gaussian pdf Is derived in this section. The inverse Gaussian pdf is defined by equation 2.39 in terras of the parameters c Q the wave velocitv, and 0, a diffusion constant (Johnson and Kotz 1970, p. 137). This form of the inverse Gaussian pdf follows naturally from the physically based parameters of the
PAGE 59
42 SaintVenant equations. Equations 2.38 and 2.40 define Cq and D, respectively, in terms of the reference discharge and cross sectional properties. However, in the ensuing analysis, the system impulse response will be determined from its moment properties which will be deduced from the changing moments of the inflow and outflow, using equations 2.56 and 2.58. To facilitate the determination of the inverse Gaussian pdf from its moment properties, an alternative form of equation 2.39 can be defined in terms of two different parameters: u, the time mean of the distribution, and X, an inverse measure of dispersion (Johnson and Kotz 1970, p. 138). Equation 2.39 can be rewritten in terms of \i and X as h(x,t)H f IG (tu,X) =/^^7^exp{X.(ty) 2 /(2.y 2 .t)} (2.63) The parameters of equations 2.39 and 2.63 are related by u = JL (2.64) c k x 2 (2.65) The mean (y IG ) f variance (Ojq), and skewness (Y iq) of the inverse Gaussian pdf are simple functions of the parameters u and X (Johnson and Kotz 1970, pp. 139140). ^IG 2 =p3 /X (2.67) IG IG
PAGE 60
43 These simple relationships facilitate the determination of state space model parameters. This is demonstrated in the following sections. Criteria for Matching Impulse Response Functions The value of a linear system analysis of the unsteady flow equations is that is allows the development of a simple model that has an impulse response function similar to the impulse response of the full SaintVenant equations. For small perturbations about a reference flow, q equation 2.63 defines a system impulse response which is the same as the impulse response for the full SaintVenant equations. In Chapter III, we will see that equation 2.63 is an excellent approximation to the SaintVenant equations for small (10 percent of Q Q ) inflow transients. The linear approximation deteriorates for large (2000 percent of Q ) inflow hydrographs. The results are also presented in Chapter III to help show the limitations of the linear theory. The strategy employed to determine the limitations of the linear theory is depicted in Figure 2.5. The primary thesis of this work is that the fundamental nonlinearities present when routing large inflow hydrographs with the SaintVenant equations can be accounted for in a model that matches the impulse response of the SaintVenant equations at all reference flow levels, Q Â• There are a number of possible criteria which could be used to compare impulse response functions for the simple and complete models. In general, functional relationships for the impulse responses cannot be obtained, so the criteria for matching impulse response functions must be able to compare the responses of the systems to impulse inputs. Some of the methods available for comparison of impulse responses are
PAGE 61
44 SRINTVENRNT EQUATIONS TRANSIENT INFLOW HYDROGRAPHS NUMERICAL SOLUTION (DWOPER) L OUTFLOW HYDROGRAPH (NUMERICAL) 1 ESTIMATION OF IMPULSE RESPONSE FUNCTION PARAMETERS \ LINEARIZED SAINTVENANT EQUATIONS 1 ANALYTICAL IMPULSE RESPONSE FUNCTION i COMPARISON OF IMPULSE RESPONSE FUNCTION 7 / UNDERSTANDING OF THE LIMITATIONS OF THE LINEAR THEORY Figure 2.5 Analysis of Limitations of Linear Theory
PAGE 62
45 1. least squares regression, 2. maximum likelihood analysis, 3. spectral analysis, and 4. matching moments. Direct derivation of parameters for the least squares criterion may not be easy. For the simplest of models, a functional relationship may be found between the parameters and the criterion to be minimized. Differentiation of that function for each parameter and solution of the set of the simultaneous equations generated is almost always difficult because of nonlinearities (Dooge 1970, pp. 170171). Statistical packages are available to perform least squares regression and to determine parameter values but, for any except the simplest models, those packages can be costly to use. Maximum likelihood parameter identification suffers from the same weaknesses. It may be difficult to determine the likelihood function for nonlinear models and the cost of using statistical packages to identify parameters involves significant cost for complex models. Spectral analysis, which involves transformations into the frequency domain, comparison of selected properties, and transformation back into the time domain, in principle, can be carried out to determine parameter values. The state space and full systems can be equated by matching pole and zero locations of the transfer functions in the transform domain. Since most hydrologic situations are not described functionally, the transformations to and from the spectral domain must be carried out numerically. Numerical inversion of the transform back to the time domain is almost always difficult (Dooge 1973, p. 31). The autocovariance properties of the impulse response function previously
PAGE 63
46 have been used to find the parameters of a state space model for a unit hydrograph (Goldstein and Larimore 1980). This approach was not pursued, but may prove useful for determining state space model parameters when the impulse response is not a singlepeaked function. The mean and variance of the impulse response functions for the linear and full systems can be compared to determine parameter values by the method of moments. If the first n moments of two impulse responses are identical, the two systems will produce identical output for polynomial input of degree n or less (Dooge 1970, pp. 170171). The method of moments is used for subsequent analysis of the state space model because it is simple and allows the convenient parameterization of the model in terms of the moments of the impulse response function. Computing the Parameters of a Gamma PDF A cascade of linear storage reservoirs is a wellstudied routing model (Zoch 1934, Nash 1959, Dooge 1973, p. 176). This model is amenable to state space formulation (SzollosiNagy 1981) and can be functionally represented by a 3parameter gamma pdf (Johnson and Kotz 1970, p. 166) k.((ta)k) n ~ 1 exp{(ta)'k} (2>69) c(t) ^ where a = the pure delay, n = the number of reservoirs, k = the time delay in each reservoir, and r() = the gamma function (Abramowitz and Stegun 1964, p. 255). As formulated in equation 2.69, n can take on any positive value. The interpretation of n as the number of reservoirs restricts it
PAGE 64
47 to positive integers, implying that T(n) could be replaced by (n1)! (the factorial function) in 2.69. Since it is useful to work with continuous parameters, the positive integer restriction on n will be relaxed temporarily. The concept of integer values for n will be reinstated for implementation of the state space model. Equations for the mean (M G ), variance (Oq), and skewness (y q) of the 3parameter gamma pdf are U G a + n/k (2.70) 2 = n/k 2 (2.71) Y = 2//k G (2.72) (Johnson and Kotz 1970, pp. 167168). Parameters of the inverse Gaussian and gamma pdfs can be related by equating their means, variances and skewnesses as defined by equations 2.66 through 2.68 and 2.70 through 2.72. When the gamma pdf parameters are found in terms of the inverse Gaussian parameters, the surprisingly simple results are a^.y (2 73) = 4_ X_ (2.74) n 9 H k 2 L (275) T M 2 As demonstrated in Figures 2.6 and 2.7, an inverse Gaussian pdf can be approximated well with a gamma pdf using equations 2.73 through 2.75. The relationships based on matching moments begin to deteriorate as n approaches 1, because the gamma pdf becomes simply a shifted exponential (Johnson and Kotz 1970, p. 166). Results for this case are shown in Figure 2.8.
PAGE 65
48 3PARAMETER GAMMA PDF DASHED LINE WITH SYMBOLS 1.0 n 4.0 k 2.0 INVERSE GAUSSIAN PDF SOLID LINE M 3.00 x 27.00 Figure 2.6 Comparison of the Inverse Gaussian and Gamma PDFs with Parameters Matched by Method of Moments for a = 1, n = 4 and k = 2.0
PAGE 66
49 3PARAMETER GAMMA PDF DASHED LINE WITH SYMBOLS a 8.0 n 8.0 k 0.5 INVERSE GAUSSIAN PDF SOLID LINE M 24.00 >> 432.00 0.0 10.0 20.0 30.0 40.0 50.0 TIME Figure 2.7 Comparison of the Inverse Gaussian and Gamma PDFs with Parameters Matched by Method of Moments for a = 8, n = 8 and k = 0.5
PAGE 67
50 3PRRflMETER GAMMA PDF DASHED LINE WITH SYMBOLS a 1.0 n 1.0 k 0.5 INVERSE GflUSSIRN PDF SOLID LINE M 3.00 x 6.75 12.0 TIME Figure 2.8 Comparison of the Inverse Gaussian and Gamma PDFs with Parameters Matched by Method of Moments for a = 1, n = 1 and k = 0.5
PAGE 68
51 Parameter values for a cascade of linear reservoirs model can be computed from values of y and X with equations 2.73 through 2.75. Values of M and X can be determined from channel properties with equations 2.64 and 2.65 if x, c and D are known, or from the mean and variance of a channel response function with (2.76) U 3 X = 11 (2.77) a 2 IG (Johnson and Kotz 1970, pp. 139140). Equations 2.76 and 2.77 can be substituted into equations 2.73 through 2.75 to provide relationships for the parameters of the gamma pdf in terms of the mean and variance of the system impulse response function (i.e., the inverse Gaussian pdf) as a T W IG M 2 n = 1. .11 (2.79) v ~ 2 IG (2.80) 3 a 2 IG In practice the mean (u_J and variance (a 2 ) would be estimated from XG Â•* J the moments of the primary model input and output hydrographs To be strictly correct in this case, the inverse Gaussian mean and variance in equations 2.78 through 2.80 should be replaced by u Z 1Q and a 2 ^, respectively, where the indicates a value estimated from data. The goal of this approach is to emulate a calibrated primary routing model (such as the numerical solution of the SaintVenant equations) with a linear state space model. The general approach is to
PAGE 69
52 1. simulate with the primary model using an inverse Gaussian pdf as the input hydrograph, 2. compute the mean and variance of the input and output hydrographs, 3. estimate the mean (jiJ and variance [o\ G ) of the system impulse response with equation 2.56 and 2.58, and 4. compute a, n and k from equations 2.78 through 2.80. When the values of a, n and k computed in step 4 above are substituted into equation 2.69, G(t) defines the impulse response of the desired linear system. A state space model with an impulse response identical to equation 2.69 is derived in the next section. A Cascade of Linear Reservoirs State Space Model The state space model developed in this section is based conceptually on a cascade of linear reservoirs. Figure 2.9 provides an overview of the steps taken to develop the linear state space model. The relationship among inflow, outflow and storage for a reservoir is described by the lumped continuity equation g!io (2 81) dt where S = the reservoir storage, I = the reservoir inflow, and = the reservoir outflow. A reservoir is linear when the relationship between storage and outflow is linear as S=<0 < 2 82) where < = a constant.
PAGE 70
53 RESERVOIR CONTINUITY EQUATION AS/At I LINEAR RESERVOIR EQUATION S K0 GENERAL STATE SPACE EQUATION Xfi E>Xt fclLai CHANNEL IMPULSE RESPONSE FUNCTION AT REFERENCE FLOW LEVEL EQUATIONS FOR LINEAR RESERVOIRS IN A CASCADE EXPRESS F AND G IN TERMS OF a, n AND k COMPUTE a,n,k AT REFERENCE FLOW LEVEL 7 STATE SPACE MODEL WITH F AND G AS FUNCTIONS OF CHANNEL" IMPULSE RESPONSE Figure 2.9 Development of the Linear State Space Model
PAGE 71
54 The state space model will operate with discrete time steps of At when implemented, so rewriting equation 2.81 for discrete time yields JiI0 (2.83) At where T = (1P )Â• I + P' I 1 2 = (lp)'O + P'O 1 2 Subscripts 1 and 2 indicate the beginning and end of a time step, respectively, and p can take on values from to 1 ; I and are weighted average values of the inflows and outflows, respectively, at the beginning and end of a time step. Substituting for T and in equation 2.83 gives 7^(lp)I +PI (lp)O p'O (2.84) At 12 12 Rewriting AS in terms of storage at the beginning and end of a time step and substituting for S with equation 2.82 produce k0 <*0 2_ L=(lp).i +pI (lp).O P'O (2.85) At 12 12 K'O + P'O 'At = Af((lp)'l + P'T ) + 2 2 1 2 (2.86) K'O (lp)'O At 1 1 or, solving for 2 (k (lp)At) 2 jrrhn r ^ 1 ** 1 !*'' 1 + mp^o > i (2 87) Equation 2.87 applies for any reservoir in a cascade. The form of equation 2.87 varies slightly for the first reservoir because the
PAGE 72
55 inflow, I t is externally specified for the first reservoir. For all subsequent reservoirs, the inflow is equal to the outflow from the previous reservoir, as I =0 ; i > 1 (288) t,i t,il where the first subscript indicates the time step and the second subscript indicates the relative position of the reservoir in the cascade. Rewriting equation 2.87 with this double subscript notation for the first reservoir yields [k (lp)At) Vi.i Tta^t + 'hJ + 1 't,i (2 89) where I t s the inflow to the first reservoir at time t and E, = k + p*At. The governing equation for all subsequent reservoirs is Vi,i=^ (H),0 t,ii +p, Vi,iJ + (2.90) [k (lo)At) 1 '"t.i where i > 1 Recall the general form of a state space model from equation 1.1 X = F X + G *U ^t+1 X,t t X,t t+1 Equations 2.89 and 2.90 must be organized into the matrix form indicated by equation 1.1. For a linear system the state transition and input coefficient matrices are constant, so equation 1.1 can be simplified to X = F'X + G'U (2.91) ^t+1 1 1+1
PAGE 73
56 Defining the state vector for a cascade of n linear reservoirs as X. = t,l t,2 t,n (2.92) and the system input vector as UP)i t + Pi t+1 Setl (2.93) F and G_ can be determined by substituting equations 2.89, 2.90, 2.92 and 2.93 into equation 2.91. Defining '< (lp).At) = T T = At T a = (ip) + p* (2.94) (2.95) (2.96) the F and G matrices for a cascade of n identical linear reservoirs are F = Tn .2 P'T^'ft n2 n1 p *T 'U y o t Â• n o pT '71 T Q (2.97)
PAGE 74
and 57 G = P'T' p n_1 .T n (2.98) Parameters k (equation 2.80) and < (equation 2.94) are related by k = 1/k (2.99) The parameter a (equation 2.69) is accounted for in the state space model by delaying the time at which the outflow from the n tn reservoir appears by a units of time. Summary A linear state space model that emulates the behavior of the full SaintVenant equations about a reference flow level has been developed. Parameters of the state space model are determined by matching the mean and variance of the system impulse responses of the primary and surrogate models. Verification of the cascade of linear reservoirs state space model by comparing the results with equation 2.39 (the impulse response of the linearized SaintVenant equations) and a study of the limitations of the linear state space model for large inflow hydrographs are the major topics of Chapter III. The work of Chapter IV, in which the important nonlinearities are analyzed and a nonlinear state space model is developed, is motivated by the problems arising when large inflow hydrographs are routed with the linear state space model.
PAGE 75
CHAPTER III LIMITATIONS OF THE LINEAR THEORY Introduction The linear theory developed in Chapter II for prismatic channels is applicable strictly only for small perturbations around a reference flow level. In this chapter, outflows predicted by the linearized SaintVenant equations and the linear state space model are compared with numerical solutions of the full SaintVenant equations for both small (10 percent of baseflow) and large (2000 percent of baseflow) inflow hydrographs. The model used as the numerical solution of the SaintVenant equations is the National Weather Service Dynamic Wave Operational (DWOPER) Model (Fread 1978). This dynamic wave routing model is based on an implicit finite difference solution of the complete onedimensional SaintVenant equations of unsteady flow. 8(A+A ) 9Q + _0 q=0 (3.1) dX d t jt+ ^__ +g .A.(i + S f ),.v^ f .B0 (3.2, where = discharge, A = cross sectional area, A = offchannel storage area, q = lateral inflow or outflow, 58
PAGE 76
59 x = distance along the channel, t = time, g = the acceleration of gravity, h = the water surface elevation, v x = velocity of lateral inflow in the xdirection, Wf = the wind term, B = channel top width, and g f = the friction slope defined as "MOM (3.3) s f 2 Â— S7T 2.2' A 'R in which n = the Manning roughness coefficient, R = the hydraulic radius, and   = the absolute value function (Fread 1978). Offchannel storage, lateral inflows and wind effects are ignored in this work. The solution technique in the DWOPER model is implicit; therefore, the time step size can be selected based on accuracy rather than numerical stability. This makes the DWOPER model very efficient in the use of computer time. Simulation of actual river reaches typically requires from 10 to 30 seconds of CPU time on an IBM 360/195 computer (Fread 1978). Verification of the Linear SaintVenant Equations The impulse response function of the linearized SaintVenant equations was derived in Chapter II as equation 2.39. (xc t) 2 h(x,t) ^ expj^ngp Â— J WD' t 3
PAGE 77
60 where and Q 2'B S ^ ; = m* v In this section, results obtained by simulating with the DWOPER model and by solving equation 2.39 are presented for 10 percent of baseflow transients on baseflows of 10000, 100000 and 200000 cfs. Inverse Gaussian shape hydrographs with a in = 3, 6, 12 and 24 hours were routed through a rectangular channel (width = 100 feet) with the DWOPER model. Transients for the same three baseflows, but for a in 6 and 12 hours were routed through a triangular channel with side slopes of 1:1. The bottom slope for both channels was 1 foot per mile. Comparisons between the numerical solution of the full SaintVenant equations and the analytical solution of equation 2.39 were made at 100 and 400 miles from the upstream end of the channels. The values of c and D in equation 2.39 were found from the channel properties at each baseflow level. Table 3.1 shows the combinations of baseflow, a* and channel geometry that produce the results for small inflow transients. It should be emphasized that we are comparing results from the numerical solution of a set of quasilinear partial differential equations with results from the analytical solution of the inverse Gaussian pdf The area under the curve defined by equation 2.39 will be the same for all values of x. The area under a hydrograph is the volume of water in the transient. Because there are no gains or losses of water in the prismatic channels studied, the volume must be equal for all x. This is guaranteed by equation 2.39.
PAGE 78
61 Table 3.1 Organization of Figures 3.1 through 3.18, Comparing the Linearized and Full SaintVenant Equations Reference Flow (1000 cfs) 10 100 200 Rectangular Channel a in (hours) 12 24 Figure 3.1 Figure 3.4 Figure 3.7 Figure 3.10 Figure 3.2 Figure 3.5 Figure 3.8 Figure 3.11 Figure 3.3 Figure 3.6 Figure 3.9 Figure 3.12 10 Reference Flow 100 (1000 cfs) 200 Triangular Channel in ( hours ) 6 I 12 Figure 3.13 Figure 3.16 Figure 3.14 Figure 3.17 Figure 3.15 Figure 3.18
PAGE 79
62 Figures 3.1, 3.2 and 3.3 show results for o in = 3 hours for baseflows of 10000, 100000 and 200000 cfs, respectively. Outflows simulated with the DWOPER model and the solution of the inverse Gaussian pdf agree closely with respect to the timing of the hydrographs for 100 and 400 miles from the channel inlet, but differ in the magnitude of the flows. In these cases the discrepancies are caused by problems with the numerical solution of the SaintVenant equations. For an abrupt wave, such as the a in = 3 hour inflow hydrograph, the distance and time step sizes for the numerical solution must be small. For the results presented in Figures 3.1 through 3.3, the distance step was 2.5 miles and the time step was 0.25 hours. With so many distance and time steps over 400 miles and 100 hours, numerical roundoff caused the DWOPER model to not maintain continuity. Therefore the flows predicted by the DWOPER model differ from those found with equation 2.39. However, the timing properties of the very abrupt inflow hydrograph are excellent. Slightly less severe inflow hydrographs, a in = 6 hours, are depicted in Figures 3.4, 3.5 and 3.6. Results from the inverse Gaussian pdf and the DWOPER model compare very well, especially for the higher flows. The numerical roundoff problems of the DWOPER model are much less evident for this inflow hydrograph. For the longer duration inflow, with o in 12 hours, presented in Figures 3.7, 3.8 and 3.9, the numerical solution of the full SaintVenant equations is matched well by the inverse Gaussian pdf. For these results, as well as those shown in Figures 3.10, 3.11 and 3.12 for a =24 hours, the analvtical and numerical solutions are almost u in indistinguishable. For these broader hydrographs the distance and time step sizes can be increased so that the numerical roundoff is no longer a problem.
PAGE 80
63 :_J :0 1Z\ Q cq Ho:: >Â•: 33 O O Â•CxJ JO CjJ!LJ! acs= ojo] i i 1 OJitu! c/rc/y OjO cncns en en; jij zd:=d: q;ql.; rt =L
PAGE 81
64
PAGE 82
65 ^ en LJ CD CD o J o 3 o OIO I I I cnco 0:0 Q_iQ_ colcn : LJItiJ COICOl Z3 ZD! Q_IQ_ Â• 0) OQ E Â•8 2 CV2 a 72
PAGE 83
66
PAGE 84
67
PAGE 85
68
PAGE 86
69
PAGE 87
70 (sjo ooo t) 3ayuHosia
PAGE 88
71 (sjo ooo n 39dyHosia
PAGE 89
72 X ^ QQ az 2: a: QJ a) CD o J a 3: a :coicoco
PAGE 90
73
PAGE 91
74 (SJO OOOI) 3QyUH0SIG
PAGE 92
75 To demonstrate the ability of the inverse Gaussian pdf to match the DWOPER model for a general prismatic channel shape, defined by equation 2.23, results for a in = 6 hours in a triangular channel are presented in Figures 3.13, 3.14 and 3.15. The delay in the arrival of the peak flow and the increased dispersion caused by the triangular channel shape can be seen by comparing, for example, Figure 3.4 with Figure 3.13. Simply by changing c Q and D of equation 2.39 to account for the new channel shape, the inverse Gaussian pdf remains an excellent approximation to the DWOPER model when the inflow transient is small. Figures 3.16, 3.17 and 3.18, for a n = 12 hours in a triangular channel, can be compared with Figures 3.7 through 3.9. Again the timing and dispersive properties of the channel are identified well by the inverse Gaussian pdf. Since the impulse response function of a linear system is a function of the system not of the inputs, cl. and D have the same values for a given baseflow for all rectangular channel results. Similarly, for the triangular channel, c Q and D are functions only of baseflow, not f in It is helpful in the analysis of, for example, Figures 3.1, 3.4, 3.7 and 3.10 (each of which is for Q = 10000 cfs in a rectangular channel), to realize that only x Q a s defined in Figure 2.4, is changing. The sets of results for the same baseflow level for either the rectangular or triangular channel shape can be interpreted as an impulse response at everincreasing distances from x = 0, the location of the upstream impulse. Figures 3.1 through 3.18 demonstrate that, for a wide range of flow levels and hydrograph durations in prismatic channels defined by equation 2.23, the impulse response function of the linearized
PAGE 93
76 io JO Xi 31 o!oi UJ.LJ: in
PAGE 94
77 x: LJ S CS o LJ CE :coicoco
PAGE 95
78 :co> O! :0 a a: 25 a: [iJ a; Â—3 Q 3Ci O J: O' 3: C2\ co;co 00 coco 00 CL.QCO CO QÂ£CC CO CO _J J Â•i" (sjo 000 n 3syyHosio
PAGE 96
79 X _, Q_ Q cc ^ en CÂ£) o LJQC Â— 1 O co>co'CO
PAGE 97
80 LJ q* 3 o L:coicoico
PAGE 98
81 3: cs o 3: o 2: :tOi
PAGE 99
82 SaintVenant equations (i.e., the inverse Gaussian pdf) produces results almost identical to those from the DWOPER model for small inflow transients. Since results from the inverse Gaussian pdf are found by analytically solving equation 2.39, the CPU time required is small (i.e., less than a second). The CPU times for the numerical solution of the SaintVenant equations on a PRIME 750 computer are presented in Table 3.2. The DWOPER model used from 1 to 3 orders of magnitude more CPU time than did the impulse response function to predict the flows shown in Figures 3.1 through 3.18. In the next section the impulse response of the linear state space model is compared with the inverse Gaussian pdf. Verification of the Linear State Space Model for Small Inflow Transients Comparison with the 3Parameter Gamma PDF The three parameters of the linear state space model specified by equations 2.91 through 2.98 are a, the pure delay; n, the number of reservoirs; and k, the time delay in each reservoir. The parameters of the 3parameter gamma pdf (equation 2.69) can be interpreted identically to the parameters of the state space model. In this section, results from the linear state space model and the 3parameter gamma pdf are compared when the same values of a, n and k are used as parameters for both the state space model and the pdf. Figure 3.19 depicts the strategy used to verify the linear state space model. The parameters c and D of the inverse Gaussian pdf can be determined for a specified prismatic channel with known slope and shape
PAGE 100
83 Table 3.2 CPU Times for the Numerical Solution of the SaintVenant Equations on a PRIME 750 Computer for Small Inflow Transients 10 Reference Flow 100 (1000 cfs) Rectangular Channel a in ( hours ) 12 24 290.15 sec 156.65 sec 49.39 sec 20.90 sec 215.23 sec 137.00 sec 40.50 sec 15.82 sec 200 I 182.50 sec 127.08 sec 39.05 sec 15.75 sec I 10 Reference Flow 100 (1000 cfs) 200 Triangular Channel o ln (hours) 6 I 12 144.56 sec 39.16 sec 141.76 sec 30.74 sec 128.87 sec 30.66 sec
PAGE 101
84 GIVEN PRISMflTIC CHANNEL APPROXIMATE THE IMPULSE RESPONSE FUNCTION WITH INVERSE GAUSSIAN PDF APPROXIMATE THE INVERSE GAUSSIAN PDF WITH 3PARAMETER GAMMA PDF SELECTED INFLOW TRANSIENTS i_Z NUMERICAL SOLUTION OF THE SAINTVENANT EQUATIONS (DWCPER) \
PAGE 102
properties. This pdf is used to approximate the impulse response of the SaintVenant equations of unsteady flow. Because the functional form of the inverse Gaussian pdf is not easily expressed in state space form, it is approximated by a 3parameter gamma pdf. Parameters of the two pdfs are related by matching their means and variances. The degree to which the 3parameter gamma pdf can approximate the inverse Gaussian pdf has been shown in Figures 2.6 through 2.8. Now we demonstrate that the linear state space model matches the gamma pdf when their parameters are equal. The state space model (like any numerical solution) cannot properly simulate given an impulse input. The impulse must be averaged over a computational time period, At. Although At can be made very small, that increases the CPU time required because it involves simulating many time steps. To avoid this problem while demonstrating that the state space model approximates well the 3parameter gamma pdf, a gamma pdf was used to generate a smooth inflow hydrograph. Gamma pdfs can be added if the parameter k is equal for all the distributions. The parameters a and n of the resultant pdf are found by summing the pure delays and the number of reservoirs, respectively, of the gamma pdfs being added. The parameters of the gamma pdf generating the inflow are defined as a r^ and k[ If the parameters of the linear state space model are chosen as a ra n m and k m > the expected model outflow can be represented as a gamma pdf with parameters a n and k if a =a. + a (3.3) o l m (3. A) n = n. + n o i ti k = k, = k o i n (3.5)
PAGE 103
86 Equations 3.3 through 3.5 were used to determine parameters for the comparisons shown in Figures 3.20 and 3.21. The same a, n and k values used in Figures 2.6 and 2.7, for the outflow hydrograph, are used in Figures 3.20 and 3.21, respectively. Values for the parameters a and ni for t he inflow hydrograph were arbitrarily chosen as and 2, respectively. The only restrictions on a i a nd n^ are and a 4 < a (36) i Â— o n < n (3.7) i o The close agreement between the linear state space model and 3parameter gamma pdf shown in Figures 3.20 and 3.21 demonstrates that the state space model is an excellent approximation of the 3parameter gamma pdf. Figures 2.6 and 2.7 show how well the inverse Gaussian and 3parameter gamma pdfs match when their parameters are related with equations 2.73 through 2.75. These results, and the excellent approximation of the linear state space model for the 3parameter gamma pdf seen in Figures 3.20 and 3.21, demonstrate that the linear state space model can emulate the inverse Gaussian pdf. Comparison with the DWOPER Model The next step in Figure 3.19 is to compare results from the state space and DWOPER models. Small (10 percent of baseflow) transients were routed with the state space and DWOPER models. The purpose here is to verify that outflow hydrographs from the linear state space model compare well with a numerical solution of the full SaintVenant equations when the inflow transients are small. To compare these two models, small transients were input to the DWOPER model. The mean and variance were computed for the inflow and outflow
PAGE 104
87 INFLOW GAMMA PDF DOTTED LINE a 0, n 2, k 2.0 OUTFLOW GflMMfl PDF DASHED LINE WITH SYMBOLS a 1, n 4, k 2.0 LINEAR STATE SPACE MODEL SOLID LINE a I, n 2, k 2.0 Figure 3.20 Comparison of the Linear State Space Model and 3Parameter Gamma PDF for a = 1, n = 4 and k = 2.0
PAGE 105
88 INFLOW GAMMA PDF DOTTED LINE a 0, n 2, k 0.5 OUTFLOW GflMMfl PDF DASHED LINE WITH SYMBOLS a 8, n 8, k 0.5 LINEAR STATE SPACE MODEL SOLID LINE a 8, n 6, k 0.5 TIME Figure 3.21 Comparison of the Linear State Space Model and 3Â— Parameter Gamma PDF for a = 8, n = 8 and k = 0.5
PAGE 106
89 hydrographs of the DWOPER model. Parameters a, n and k of the 3parameter gamma pdf (identical to the parameters of the linear state space model) were computed from these two moments using equations 2.78 through 2.80. These three parameters, with equations 2.91 through 2.98, define the linear state space model. Varying the value of p, the time weighting factor, from to 1 had no impact on the results in the cases tested. The parameter p was set to 1 for the results presented to assure stability of the linear equations (Fread 1974, pp. 710). For purposes of this comparison, selected channel geometry and inflow combinations from the previous section were repeated. Table 3.3 shows pairs of figures with identical inflow and channel properties. Now, reimpose the integer restriction on parameter n. In the state space model, n represents the number of reservoirs in the cascade and must, therefore, be a positive whole number. From equations 2.70 and 2.71, the mean and variance of the 3parameter gamma pdf (and also of the impulse response function of the state space model) are U = a + n/k and a 2 = n/k 2 The exact values of a, n and k used in the state space model are presented in Table 3.4. These values were computed from the moments of the inflow and outflow hydrographs of the DWOPER model using equations 2.78 through 2.80. However, n must be an integer. The values of n listed in Table 3.4 were rounded to the nearest integer for use in the state space model. The slight effect of this parameter change can be seen by comparing Figure 3.22 with Figure 3.4. The ideal value of n is 3.702,
PAGE 107
90 Table 3.3 Pairs of Figures with Identical Input Data Comparing the Linearized SaintVenant Equations and the Linear State Space Model with the Full SaintVenant Equations Figure Comparing Figure Comparing Linear Linearized and Full State Space Model and Full SaintVenant Equations SaintVe nant Equations 3.4 3.22 3.2 3.23 3.8 3.24 3.12 3.25 3.14 3.26 3.18 3.27
PAGE 108
91 Table 3.4 Values of Parameters a, n and k for the State Space Model Rectangular Channel, L = 100 miles
PAGE 109
92 en a: cs 9Â§ CS
PAGE 110
93 but this must be rounded to 4. Equations 2.70 and 2.71 indicate that this change should increase the lag and the dispersion introduced by the state space model. This is exactly what occurs. Whereas in Figure 3.4, the timing of the peak for the impulse response was identical with the DWOPER model and the peak flow slightly high at 100 miles, the peak flow from the state space model is later and lower in Figure 3.22. The rounding of parameter n can sometimes seem beneficial, as Figure 3.23 shows. The ideal value of n is rounded from 0.943 to 1. This will increase the lag and the variance of the state space model outflow hydrograph. By slightly delaying the rising and the falling limbs of the hydrograph, the state space model flow agrees better visually with the DWOPER model than did the inverse Gaussian pdf in Figure 3.2. A more expected result can be seen by comparing Figures 3.24 and 3.8. In this case, the increased variance has caused the peak flow to decrease because the volume of water in the hydrograph is fixed, and the falling limb of the hydrograph for the state space model has been delayed. For a baseflow of 200000 cfs in a rectangular channel, n would ideally be equal to 0.625. The change in the outflow hydrograph from the impulse response function of Figure 3.12 to the state space model outflow shown in Figure 3.25 is caused by the substantial correction of parameter n to an integer value. Since the maximum time shown on Figure 3.25 is larger than for the other figures, horizontal differences between curves are correspondingly larger. A case of n being rounded down to the next lower integer can be seen in Figure 3.26 for a triangular channel. In this case the lag and
PAGE 111
94 en z a: J e J a o OiC z:io_ ICO O (C 3H a!cn i 05 a i o Â§ j> Q CO % I Â£ I "to os 9 o o_ c * o i i < cm a ^ CM ^ r { s ........................................ 8 5 (sjo ooo n sMOUino i3aow qnu moijni
PAGE 112
95 1=
PAGE 113
96 X
PAGE 114
97 X: o_: Si 9S\ UJ CD i !_. ! J s uu Q O l0_ uu a. ti o cc 2H olto I GO O I CO DC to t 6 CC J Â• a. c Â• o ? i oo a ^ CO o R S W CO O UJ R CD GO .SI 3 3 tsjo ooo t) sMOHJino i30ow onu moijni
PAGE 115
98 variance of the state space model decrease because n is less than the ideal. Comparison of Figure 3.26 with Figure 3.13 demonstrates that the state space model outflow is distinctly earlier and less dispersed than results from the impulse response function. This causes the peak flow for the state space model in Figure 3.26 to be higher than for the DWOPER model. State space model results for a baseflow of 200000 cfs in a triangular channel are presented in Figure 3.27. The effects of increasing n from 0.874 to 1 can clearly be seen by comparing the outflow hydrograph with Figure 3.18. The peak flow of the state space model in Figure 3.27 is lower because the outflow hydrograph is more dispersed. The increased lag can be seen in the hydrograph's falling limb. Results from the linear state space model compare very well with those from the DWOPER model. Small differences in results can be attributed to the interger restriction on parameter n. Differences in CPU times required to simulate with the two models are large. The CPU times on a PRIME 750 computer for the six figures discussed above are presented in Table 3.5. The DWOPER model typically took from 2 to 3 orders of magnitude more CPU time than the linear state space model. The Importance of Nonlinearities The previous results show that the linear theory developed in Chapter II works very well when the phenomena being modelled are linear. However, in many (if not most) of the surface runoff or channel routing problems encountered, the input fluctuations are large, and nonlinear behavior is observed. In this section the linear state space
PAGE 116
99 Â— Â— rÂ—
PAGE 117
100 Table 3.5 CPU Times for the Linear State Space and DWOPER Models on a PRIME 750 Computer for Results of Figures 3.22 through 3.27 Linear State
PAGE 118
101 model is compared with the DWOPER model for large (2000 percent of baseflow) inflow transients. Since the linear model is only strictly applicable for a small flow range around the reference discharge, we expect the results of this section to deterioriate significantly from those presented above. If this were not the case, Chapter IV would be titled "Conclusions" instead of "A Nonlinear State Space Model." The cases presented are for inflow hydrographs of 200000 cfs on a baseflow of 10000 cfs. Results from the linear state space and DWOPER models are compared at 100 and 400 miles from the channel inlet. The 16 cases presented are summarized in Table 3.6. As shown in Figures 3.28 and 3.29, a linear model is not adequate to simulate the physical processes that occur when a large, very abrupt hydrograph (o in = 3 hours) is input. Since results from the linear state space model are sensitive to the reference flow level at which the parameters a, n and k are selected, Q was chosen as 100000 cfs for all large transient figures unless otherwise stated. This value was chosen at slightly less than the average of the base and peak flow levels, to assure that values for the parameters of the state space model are as often too high as too low. The results of Figure 3.28 indicate that, even at 100 miles from the channel inlet, the state space model does not disperse the inflow transient as much as the DWOPER model. At 400 miles, shown in Figure 3.29, outflow from the state space model is clearly late and not as dispersed as from the DWOPER model. Figure 3.30 shows that the timing and dispersion properties of the state space model for o ln = 6 hours at L = 100 miles are very close to those of the DWOPER model; however, the rising and falling limbs of the state space model outflow hydrograph have the wrong slopes. This is a
PAGE 119
102 Table 3.6 Organization of Figures 3.28 through 3.43, Comparing the Linear State Space Model with the Full SaintVenant Equations for Large Inflow Transients Rectangular Channel 'in 'in = 3 hours 'in 'in 12 hours L = 100 miles Figure 3.28 Figure 3.30 Figure 3.32 Figure 3.34 L = 400 miles Figure 3.29 Figure 3.31 Figure 3.33 Figure 3.35 Triangular Channel a, 6 hours u in "i n = 12 hours L = 100 miles Figure 3.36 Figure 3.38 L = 400 miles Figure 3.37 Figure 3.39 DoublePeaked Hydrograph = 100000 cfs Qq = 10000 cfs Q = 200000 cfs L = 100 miles Figure 3.40 > = 400 miles Figure 3.41 Figure 3.42 Figure 3.43
PAGE 120
103 cc a: 9Â§ 14
PAGE 121
104 CC a: q S i !_j J S uu Q O OiCC z:io. LJ Id OJli O (C 3H oicn I en o I o Â§ 4) Qfi Â•it:::::: Â— 1 i tftt.t T ............... R 8 n 8 R 8 N 8 R (SJO 0001) SM01JM0 13Q0W QNy M01JNI
PAGE 122
105 CD T
PAGE 123
106 typical problem with linear models trying to simulate nonlinear behavior. The linear results are closer to a Gaussian pdf than to the truely nonlinear shape of the DWOPER model outflow hydrograph. This will be a recurring problem for the linear state space model results presented in the remainder of this chapter. Results for a n = 6 hours at L 400 miles, shown in Figure 3.31, demonstrate that it was merely good fortune, not a properly formulated model, that let the linear state space model match the DWOPER model outflow in Figure 3.30. At 400 miles, the timing and dispersive properties of the state space model are clearly different than the nonlinear results from the DWOPER model. With the broader inflow hydrograph (o in = 12 hours) shown in Figure 3.32, the behavior of the state space model improves greatly. Since the flow varies more gradually over time, the nonlinear terms of the SaintVenant equations are relatively less important than for the more rapidly varying inflows shown earlier. The major discrepancies between the state space and DWOPER models are the slight increase in lag and dispersion of the inflow hydrograph shown by the state space model which is caused by rounding parameter n up to 1, and the symmetrical shape of the state space model's outflow hydrograph, which is typical of linear models. The Gaussian shape of the state space model outflow is even more evident in Figure 3.33, for a n = 12 hours and L = 400 miles. For this slowlv varying inflow, the delay and dispersion introduced by the state space model are approximately equal to those seen in the DWOPER model results, but the linear model is not able to reproduce the slopes of the steeply rising and slowly falling outflow hydrograph limbs which are characteristic of nonlinear routing models.
PAGE 124
107 a:
PAGE 125
108 CSJO 0001) SMOUinO 13Q0W QNU MOHJNI
PAGE 126
109 X
PAGE 127
110 As the dispersion of the inflow hydrograph increases, these effects become even more pronounced. In Figures 3.34 and 3.35, for L = 100 and 400 miles, respectively, the basic timing and attenuation properties of the state space and DWOPER models are in close agreement, but the symmetrical shape, typcial of the linear model, limits its ability to match the nonlinear outflow hydrograph. For a triangular channel, as shown in Figure 3.36, the major cause of differences between the state space and DWOPER models is the rounding of parameter n down to 1. As we will see in Chapter IV, the correct value for m, the kinematic wave parameter, is 4/3 for a triangular channel. This is close to m = 1, the value implied by a linear model. Since m should be equal to 5/3 for a rectangular channel, we would expect a linear model to be better for a triangular channel than for a rectangular one. This is indeed the case, as can be seen by comparing Figure 3.37, for a^ n = 6 hours and L = 400 miles in a triangular channel, with the corresponding Figure 3.31 in a rectangular channel. The shape of the state space model outflow in Figure 3.37 is closer to the DWOPER model than in Figure 3.31. This is not because the state space model is any less linear in Figure 3.37, but because flow in a triangular channel, with m = 4/3, is more linear than in a rectangular channel, where m = 5/3. This effect is even more evident in Figures 3.38 and 3.39 for an inflow hydrograph with a ln = 12 hours at L = 100 and 400 miles, respectively. The state space model outflow hydrograph shapes, especially the falling limbs, are much closer to the nonlinear behavior simulated with the DWOPER model than in the corresponding Figures 3.32 and 3.33, for a rectangular channel.
PAGE 128
Ill (sjo oooi) SMOUino i3Qow qnu mum
PAGE 129
112 L k CO ZD O 3: Â•5 =C re iÂ— i Ls is ha 8 ~"s "^ S 8 8 ~TT 8 8 ol w ** "* Â•*Â• *"Â• (sjo ooo n sMOumo naoow gnu moijni u o o
PAGE 130
113 jo.
PAGE 131
114 GO 00= s CSC o Q J 2 LJ LJ Q O OiCC niann uj;lj a. h O (C a! en i i rÂ§ K 8  en a: Â•si o ^T 78 K Â•8 8 V T" K "T 8 8 Â§ S S 8 Â§ CSJO 0001) SHOTJlflO 1300W ONO MOUNI
PAGE 132
115 CSJ3 0001) SMOIJinO 1300W om motjni
PAGE 133
116 m
PAGE 134
117 The next test comparing the linear state space and DWOPER models was for a doublepeaked inflow hydrograph. The inflow hydrograph consisted of a 100000 cfs transient followed by a 200000 cfs transient on a baseflow of 10000 cfs. At 100 miles from the inlet of a rectangular channel, the linear state space model still approximates the DWOPER model results well. However, as shown in Figure 3.40, the rising and falling linbs of the linear model's hydrograph are already different in slope than the DWOPER model. In Figure 3.41 we can see that, by 400 miles from the channel inlet, the state space model outflow has departed significantly from the DWOPER model. The outflow hydrograph from the linear model is becoming more symmetrical, while the rising limb of the DWOPER model outflow is steepening and the falling limb is becoming flatter. The impact of the reference flow level on the linear state space model results is shown in Figures 3.42 and 3.43. The parameters a, n and k were selected at a reference flow of 10000 cfs in Figure 3.42. The resultant outflow hydrograph is very late and much more dispersed than the flows shown in Figure 3.41. When parameters are chosen at q = 200000 cfs, as in Figure 3.43, the state space model lags and disperses the inflow hydrograph significantly less than when Q = 100000 cfs. The reference flow level has a striking effect on the linear state space model behavior. Summary Results presented in this chapter have demonstrated 1. the behavior of the impulse response function for the linearized SaintVenant equations (equation 2.39) is almost
PAGE 135
118 a
PAGE 136
119 31
PAGE 137
120
PAGE 138
121 (SJO OOQI) SMOUinO H3Q0W GNU M01JNI
PAGE 139
122 identical with a numerical solution of the full SaintVenant equations for small (10 percent of baseflow) inflow transients 2. the linear state space model developed in Chapter II is an excellent surrogate for a model based on the full SaintVenant equations when small inflow transients are routed, and 3. important nonlinearities that are not adequately modelled by the linear theory developed in Chapter II when large (2000 percent of baseflow) inflow hydrographs are routed with the full SaintVenant equations. In Chapter IV nonlinearities important when routing large inflow hydrographs are analyzed. A nonlinear state space model emulate the behavior of the full SaintVenant equations over large flow ranges is developed. The ability of this nonlinear state space model to act as a surrogate for the DWOPER model is the topic of Chapter V.
PAGE 140
CHAPTER IV A NONLINEAR STATE SPACE ROUTING MODEL Analysis of the Important Nonlinearltles Chapter III presented limitations of the linear theory. The inability of the linear theory (specifically the impulse response function of the linearized SaintVenant equations) to emulate the behavior of the full nonlinear SaintVenant equations leads us to the current chapter. In this chapter we will 1. analyze the nonlinear characteristics of the primary model, 2. explore the properties of a nonlinear storage reservoir, and 3. develop a fully nonlinear state space routing model which emulates the mathematical behavior of the SaintVenant equations for any inflow transient. The approach to the analysis of the nonlinearities inherent in the SaintVenant equations is outlined in Figure 4.1. The following section begins with the full SaintVenant equations and, in the spirit of the new perspective of routing presented in Figure 1.1, predicts the theoretical variation of the impulse response properties with a linear analysis. The expectations of linear theory are then confirmed by numerically solving the fully nonlinear SaintVenant equations and observing the changes in the impulse response function. Results from this analysis lend insight into a state space model structure which will emulate the nonlinear behavior of the primary model. 123
PAGE 141
124 FULL NONLINEAR SRINTVENRNT EQURTIONS PREDICT HOW MOMENTS OF IMPULSE RESPONSE FUNCTION VARY WITH FLOW LEVEL 1ENTS Zi NUMERICAL SOLUTION (DWOPER) Z OBSERVE HOW MOMENTS OF IMPULSE RESPONSE FUNCTION VARY WITH FLOW LEVEL LZ MOMENTS OF IMPULSE RESPONSE VARY AS ALGEBRAIC POWER OF FLOW LEVEL I CONSTRUCT NONLINEAR STATE SPACE MODEL USING CASCADE OF NONLINEAR RESERVOIRS AND VARIABLE LAG L COMPARE RESPONSE OF NONLINEAR STATE SPACE MODEL WITH RESULTS FROM PRIMARY MODEL Figure 4.1 Approach to the Nonlinear Analysis
PAGE 142
125 Nonlinear Characteristics of the Primary Model Introduction This section analyzes the nonlinear characteristics of the SaintVenant equations, the exemplar primary model. The approach for that analysis is shown in Figure 4.2. The variations with flow of the linearized SaintVenant equations' impulse response function are derived and the mean and variance properties of the impulse response function are analyzed. These moments will be needed to specify the parameters of the nonlinear state space model developed later. The variations with flow of the impulse response for the full SaintVenant equations are determined by simulating with incremental transients as inflow, and observing the outflow. The results predicted by the linear analysis of the SaintVenant equations agree with those obtained by simulating with the DWOPER model. Theoretical Variations with Flow Level of the Impulse Response Function for the Linearized SaintVenant Equations for Incremental Inflow Transients Beginning with the SaintVenant equations of unsteady flow (equations 2.20 and 2.21), the impulse response function for the linearized SaintVenant equations was obtained in Chapter II as equation 2.39 (xc t) 2 4'TTD't 3
PAGE 143
126 USE ANALYTICAL IMPULSE RESPONSE FUNCTION TO PREDICT VARIATION OF IMPULSE RESPONSE FUNCTION WITH FLOW LEVEL TRANSIENT INFLOW HYDROGRAPHS I NUMERICAL SOLUTION (DWOPER) i DETERMINE VARIATION OF IMPULSE RESPONSE FUNCTION PROPERTIES WITH FLOW LEVEL V7 COMPARE RESULTS I CONCLUSIONS: M wQ [lrl)/i] n [(3r2m)/2m] 9 *Â•* U Figure 4.2 Analysis of the Nonlinear Characteristics of the Primary Model
PAGE 144
127 This equation has the form of an inverse Gaussian pdf (Johnson and Kotz 1970, p. 137). The mean, t^, and variance, a 2 of the inverse Gaussian pdf were presented in equations 2.41 and 2.42 as and 'if a 2 = 2 Â— fiV S c 3 The derivation of the changes with flow level, 0, in the mean and variance of the linearized SaintVenant equations' impulse response function begins by substituting for c in equation 2.41 with equation 2.38 to give t. *(4.1) x. ra* v From equation 2.22 so t. =Â•(42) I m Equation 2.30 expresses the friction slope as 2 S =_^ S f 2 2m n a 'A u Low Since this analysis is for small transients about a reference flc level, we make the standard kinematic assumption of steady uniform flow S c = S (4.3) f (Henderson 1966, pp. 287,367). Substituting equation 4.3 into equation 2.30 and rearranging gives
PAGE 145
(i) 1/m Q 1/m Replacing A in equation 4.2 with equation 4.4 produces 1/m s x, f M/B,Q A m ^^ 128 (4.4) (4.5) 1m h = v Q where x f U 1/m Recognizing that 8. is constant for a reference flow 1m t t ~ (4.6) (4.7) where ~ symbolizes "is proportional to." Equation 4.7 states that the first moment of the impulse response function varies as a simple algebraic power of the flow level. The value of the exponent depends upon the kinematic wave parameter, m. The variation of the second moment can be found by again substituting for c with equation 2.38. This time substitution into equation 2.42 gives ,2 Â„ Qx B S Â• m 3 Â• v 3 1H>' (4.8) where f b q and c Q are now replaced by Q, B and c because we are solving for the variation of a 2 with flow. Substituting for v with equation 2.22 produces B'S m 3 Q 3 (4.9)
PAGE 146
129 By replacing B with equation 2.23 we get a 2 = 9^. ,(*!).( i*2) (4.10) B A r S Â•m 3 3 And, finally, substituting for A with equation 4.4 and rearranging yields 3r Â£2t2 Â„ I Â• S Â• m 3 (iV)Ci) m '^Â— (4.11) By simplifying the expression for Q we obtain 3r2m a 2 = Q m (4.12) 2 where m e = l_(i* 2 )a r 3 2 B Â• S Â• ra 3 The standard deviation of the impulse response is 3r2m a =/r /bq 2m (4.13) 3r2m a ~ Q 2m (4.14) Again, the moment of the impulse response varies as a simple algebraic power of and the exponent depends on m. For the variance, however, the power is also a function of the channel shape parameter, r. For a given friction resistance formula (i.e., Manning or Chezy) m and r are related by equation 2.28. The coefficient of variation of the linearized SaintVenant equations' impulse response can be found from equation 2.43
PAGE 147
130 a Substituting equations 4.6 and 4.13 into equation 2.43 gives 3r2m /3 2m %=^ Q T^ (4 15) 1 Q~ or simplifying 1r c = 6 Q 2m (4.16) v 3 where /3~ Q 2 3 3 X 1 or finally 1r c ~ 2m (4.17) v Equations 4.7, 4.14 and 4.17 express the theoretical variations with flow level of the mean, standard deviation and coefficient of variation of the impulse response function of the SaintVenant equations. In the following section some numerical results are presented for comparison with the variations predicted by linear theory. Simulated Variations with Flow Level of the Impulse Response Function for the Full SaintVenant Equations for Incremental Inflow Transients In general, the operational procedure for determining the moments of the channel response will be to simulate small transients with the primary model at various baseflow levels and observe the results. In
PAGE 148
131 the current application, however, with the SaintVenant equations as the primary model for a prismatic channel, we have the opportunity to compare the expected theoretical results obtained from equations 4.7, 4.14 and 4.17, with numerical computation of the moments from the fully nonlinear SaintVenant equations model. The NWS DWOPER model (Fread 1978) was used to numerically solve the full equations of unsteady flow. The moments of the impulse response function of the SaintVenant equations were found by driving the DWOPER model with small (10 percent of baseflow) inflow transients for several baseflow levels, and observing the outflows. The means and variances of the inflow, x, and outflow, y, hydrographs were then computed. The moments of the impulse response function, h, were found with equations 2.53 and 2.54 U'(h) = u'(y) U'(x) and U (h) = U (y) U (x) 2 2 2 The coefficient of variation of the impulse response is /U (h) c = 2 (4.18) C v U'(h) 1 From equations 4.7, 4.14 and 4.17, we see that the theoretical variations with flow of the impulse response moments depend on m and r. The theoretical values of m can be derived for any prismatic channel using the Manning friction formula. Once m is known, the value of r can be found from equation 2.28. For the two channel shapes studied in this work, the values of m and r are:
PAGE 149
132 shape m r rectangular 5/3 triangular 4/3 1/2 Figures 4.3 and 4.4 depict the lag time introduced by the rectangular and triangular channels, respectively. The geometry of the channels is described in Chapter III. The lags presented are for distances of 25, 50, 100, 200 and 400 miles from the channel inlet (Xq in the notation of Chapter III). Incremental transients on baseflows of 10000, 100000 and 200000 cfs were routed with the DWOPER model. The results presented are for inflow transients with o" in = 12 hours. Results for a.. 3, 6 and 24 hours look identical to those presented. This must be, because the impulse response for a linear system is not dependent on the properties of the input, but is strictly a function of the system. The theoretically predicted and numerically simulated values match almost exactly. The slope of the lines on Figures 4.3 and 4.4 can be found by rearranging equation 4.7 as slope, Â™ (4.19) lag 1m For the rectangular channel presented in Figure 4.3, with m = 5/3, the slopei = 2.5. On Figure 4.4, the triangular channel with m = 4/3, the slope^ = 4. Observed lag time versus flow data which corroborate the theoretical expectations of equation 4.19 are presented in Chapter VII. Figures 4.5 and 4.6 present the theoretical and simulated standard deviations of the impulse response function versus flow level. Again, the results compare very well. The slopes of the lines on Figures 4.5 and 4.6 can be found by substituting (m = 5/3, r = 0) and (m = 4/3, r = 1/2), respectively, into
PAGE 150
133 az arco UJOCxJ zuz _]{Â• _j i i i i i to t~> m^Z * CM T 1 2=^^ I I I I I cc iÂ— co O 5Â— Â•Â— i aoa_ o < X o en !jj IÂ£
PAGE 151
134 CO"cccccn OJOUJ z u 2: v1I11 Â—1 OJ CO "* LT> _Jt_J I I I I CUE' UD 0_O_ Oh QO uinoooo _j oj m o o o i a a b ^O < + x o O en UJ L 5 I I I I r 1 I I 4 o 3 ISJO) 39yUH0SIQ 30N3;d333y
PAGE 152
135 rS cccccn uoy zuz _jt_j i i i i i en tÂ— _]0_JCMLnoao .Â— a: z r^o3* cti ~ i i I i i Q_ Ql. Z OHm aoa_ o < + x o en o x CO C\2 cc cc:
PAGE 153
136 o:>.~ cncccn uou .xiÂ— > *< oz m v m 3 eÂ— Â• Â— 3 uEuinoooo JQJtMl^OOO SJO) 39^UH0SIQ 30N3y3J3^
PAGE 154
2m ,lope a = T^i" 137 (4.20) Similarly, Figures 4.7 and 4.8 depict the theoretical and simulated coefficients of variation of the impulse response function. The slopes of these lines can be found as slope 2m 1r (4.21) Before leaving the topic of the simulated variance with flow level of the impulse response function, a valuable relationship can be found by solving equation 4.19 for m as slope lag 1 + slope.. lag (4.22) This expression, defining the kinematic wave parameter, ra, in terms of the channel lag versus flow properties, will be useful for determining one of the parameters of the nonlinear state space model which is developed in a subsequent section of this chapter. Properties of the Nonlinear Storage Reservoir as a Building Block for the State Space Model Just as the linear reservoir relationship, equation 2.82, is the cornerstone of the linear state space model, the nonlinear reservoir equation is the foundation in the development of the nonlinear state space routing model. The relationship between storage and outflow for a nonlinear reservoir is S = <Â• where S = the reservoir storage, (4.23)
PAGE 155
138 LINEAR THEORY LINES) 1
PAGE 156
139 LINEAR THEORY (LINES) 1
PAGE 157
140 = the reservoir outflow, and m and < are constants. If m in equation 4.23 is taken to be the kinematic wave parameter, the nonlinear storage reservoir has the same lag versus flow properties as the linearized SaintVenant equations. Equation 4.23 can be rewritten as s = < 1/m .o 1/m (4.24) For a reservoir of length L and cross sectional area A S = LA (4.25) Combining equations 4.24 and 4.25 gives 1 1/m 1/m ,. ... A = jk 0 (4.26) The time lag through a reservoir is t. =1 (4.27) Jc v where v is the average velocity through the reservoir. But from equation 2.22 V A Substituting equation 2.22 into equation 4.27 produces h = L 4 (4 28) And substituting equation 4.26 into equation 4.28 yields
PAGE 158
141 1m t = < 1/m m (4.30) which is identical in form to equation 4.6 for the time lag of the linearized SaintVenant equations. The thesis presented in this work is that a state space model in the form of a cascade of nonlinear storage reservoirs can emulate the nonlinear behavior of the SaintVenant equations. A fully nonlinear state space model that maintains the lag versus flow relationship of equations 4.6 and 4.30 is developed in the next section. Development of the Nonlinear Model Introduction A fully nonlinear state space routing model which emulates the SaintVenant equations for downstream waves in prismatic channels is developed in this section. This nonlinear state space model is based upon 1. the linear state space model developed in Chapter II, which has the same Impulse response function as the linearized SaintVenant equations for a prismatic channel (i.e., the Inverse Gaussian pdf); and 2. the functional relationships for the variations with flow of the mean and variance of the inverse Gaussian pdf which were derived earlier in this chapter.
PAGE 159
142 Figure 4.9 shows the approach taken in developing the nonlinear state space model. The general structure of the nonlinear state space model is a cascade of nonlinear reservoirs as shown in Figure 4.10. The nonlinear properties are implemented by varying, with flow level, the parameters (a, n and k) of the linear state space model developed in Chapter II. Instead of a model specified by three parameters, as In Chapter II, the nonlinear state space model is specified by three functions. The values of a, n and k vary with the level of flow. Although each of the n reservoirs in the cascade is governed by identical functions of a versus flow and k versus flow, in general the values of ai and k.< will be different for each i th reservoir because, except under steady flow conditions, the flow through each reservoir will be different. From three functional relationships, 2n+l parameters are determined which specify the F^ t an d _Gj( t matrices in equation 1.1 at each time step. The algorithm to perform this task efficiently is presented later in this chapter, but first the functions of a, n and k versus flow are derived. Equations 4.7 and 4.14 express the variations with flow of the lag and standard deviation of the impulse response function for the linearized SaintVenant equations. We want the variations with flow of the lag and variance properties of the state space model to match the variations of the fully nonlinear SaintVenant equations. The parameters of the state space model are related to the mean and variance of the impulse response function by equations 2.78 through 2.80
PAGE 160
143 ASSUME STRUCTURE OF STRTE SPRCE MODEL EQUIVALENT TO fl CASCADE OF NONLINEAR RESERVOIRS PLUS A PURE DELAY THAT VARIES WITH FLOW LEVEL DERIVE ANALYTICAL EXPRESSIONS FOR MOMENTS OF THE IMPULSE RESPONSE FUNCTION OF THE STRTE SPACE MODEL AS A FUNCTION OF FLOW LEVEL i EXPRESS PARAMETERS OF 3PARRMETER GAMMA PDF AS A FUNCTION OF MOMENTS AND FLOW LEVEL i EXPRESS STATE TRANSITION AND INPUT COEFFICIENT MATRICES AS FUNCTIONS OF FLOW LEVEL Figure 4.9 Approach to the Development of the Nonlinear State Space Model
PAGE 161
L44 UPSTREAM INFLOW PURE TIME DELAY DOWNSTREAM OUTFLOW Figure 4.10 A Cascade of Nonlinear Reservoirs Governed by S 111 = k*0
PAGE 162
145 1 a = T*
PAGE 163
146 Equations 4.31, 4.33 and 4.35 describe the variations with flow of parameters a, n and k, respectively, of the state space model. Operationally, the value of m for the channel can be found by simulation. Using the primary model, the lags for several baseflows with small input transients can be found. The slope of the line obtained when lag is plotted as a function of baseflow on logarithmic axes determines m with equation 4.22. The lag and variance produced by the channel must be found for any one flow level by simulating with the primary model. With this as a reference point in lagQ space, ra computed as above, and the functions 4.31, 4.33 and 4.35 for the model parameters, all the information needed to construct the nonlinear state space model is at hand. The steps required to develop the nonlinear state space model are summarized in Figure 4.11. Two points of a general nature should be stressed at this time 1. The state space model developed here is fully nonlinear. Its theoretical basis draws heavily from linear system analysis, but only to explain better the nature of the important nonlinearities; and 2. There is not a reference discharge for this model. The functions described by equations 4.31, 4.33 and 4.35 are valid for all flows where m is constant. A point in lagQ space to fix these functions could be found at any flow level. The following sections describe the final sucessful nonlinear state space model and the path taken to find the final model. The first is needed to confirm the perspective of this work. The second is
PAGE 164
147 RESERVOIR CONTINUITY EQUATION AS/Ai I NONLINEAR RESERVOIR EQUATION S c0 GENERAL STATE SPACE EQUATION CHANNEL IMPULSE RESPONSE FUNCTION AT SEVERAL REFERENCE FLOW LEVELS DETERMINE FLOW RANGES IN WHICH m IS CONSTANT AND COMPUTE a n AND k AT A REFERENCE FLOW IN EACH RANGE EQUATIONS FOR NONLINEAR RESERVOIRS IN A CASCADE i NONLINEAR STATE SPACE MODEL WITH F X1 AND G xt AS FUNCTIONS
PAGE 165
148 included because it may lend insight into the alternatives available even after the theory is developed. Obviously, both the theoretician and the model builder are vital if a useful product is to be developed. The Final Configuration of the Nonlinear Model The final form of the nonlinear state space routing model is based conceptually on a cascade of nonlinear reservoirs. As with the linear state space model developed in Chapter II, equation 2.81 describes the relationships among inflow, outflow and storage for a reservoir dT= I " where S = the reservoir storage, I = the reservoir inflow, and = the reservoir outflow. Now, however, the relationship between storage and outflow is nonlinear, as expressed in equation 4.23. The nonlinearites implicit in equation 4.23 are modelled by letting k vary with flow for the i th reservoir as S = K Â•() (4.36) t,i t,i t,i where the subscript t,i indicates the value of < for the flow level through the i th reservoir during the t th time step. Equation 2.99 defines parameter < as the inverse of parameter k. The variation of k with flow is defined by equation 4.35. Equation 2.84 is still valid for nonlinear reservoirs
PAGE 166
149 Â£! (ip)i +pi (ip)o po At i 2 1 2 Rewriting AS in terms of storage at the beginning and end of a time step for the i tn reservoir yields S . S t4_ t 1 diH ttH+ p.o ttltl 1 (4.37) a p) t,ip *Vi,i Substituting for S. }i with equation 4.36 produces t+l,i t+1,1 t,i t,i = At (1 ^ )  t ,ii + p Vi,ii(1 ^ ) *t,ip *Vi,i (4.38) Solving for O t+ i s i [< (lp)At) o (iP)^ Q + l t,i ^_. + 6+1,1 C^,, + PAtj t 1 1 (< t+1 i+ P'At) (4.39) 7 ~T t+i,ii (ic t+11 + P'At) The value < t+1 i in equation 4.39 depends upon Ot+1,1To avoid the need for an iterative solution let k < 4 (4.40) t+l,i t,i Using the approximation of equation 4.40 and defining E = k + P'At (4.41) t,i t,i
PAGE 167
150 and (* tfl (lP>At] 1 t>x (4.42) t 1 5 t,i equation 4.39 becomes (lp).At +^0 ,,, (4.43) t+i.i iÂ„ t,ii t,i t,i r~7 t+i f ii t 1 C 1 For the first reservoir in the cascade, 0. j_j i s t he specified system inflow, so equation 4.43 reduces to Af((lp)I + P'I t+1 ) , * Â•() + ElL(4.44) t+1,1 t,l t,l 4 t t By defining Vi = (^>i t + P'Vi (4 45) the governing equation for the first reservoir in the cascade becomes = 0 + At 'u (4.46) t+i,i t,i t,i ftt t+i t 1 The governing equations for subsequent reservoirs can be found by recursively substituting for t+1>1 _! in equation 4.43. Define T = *L_ (4.47) t 1 ^7 and fl = (1p) + pt (4.48) t i t ,i The value of the element for the j tn row and m tn column of the j^ t matrix in equation 1.1 is
PAGE 168
151 zero for m > j (4.49) for m = j (4.50) t ,tn fl T for m = 11 (4.51) t,m t,j a .p^"Â™" 1 n T t for in < j1 (4.52) **" Wl t' 1 where II is the product operator. The ith j cri row of the G^ fc matrix i j J" 1 n T (4.53) 1=1 t,i Equations 4.49 through 4.53 complete the information needed to construct a nonlinear state space routing model. The steps that have led to this point can he summarized as 1. development of a linear state space routing model, having constant parameters a, n and k, with an impulse response identical to the linearized SaintVenant equations, 2. determination of the variations with flow level of the mean and variance of the linear state space model impulse response function, 3. relation of the variations with flow of the impulse response function's moments to variations of the parameters a, n and k, and 4. derivation of the state transition and input coefficient matrices as functions of flow level. The theory is now well defined: let the model parameters vary with flow and simulate using equation 1.1. There are, however, many
PAGE 169
152 ways to implement the theory, and many possible algorithms for simulation. The general approach in developing the nonlinear state space model is to retain as much of the theory as possible while keeping the model algorithm simple. The final algorithm, presented here, incorporates all the theory developed in this dissertation. The preliminary algorithms are described in the following section. As more of the theory was introduced into the approaches to state space modelling of unsteady flow, the results compared more favorably with the solution of the full SaintVenant equations. The final algorithm emulates well the solution of the complete equations of unsteady flow over a wide range of channel and wave conditions. The algorithm for simulating with the nonlinear state space model is: 1. Compute the average of the outflows from each of the n reservoirs at time = t, and the current system inflow Q (u + I J/(n+l) (4.54) t+l t+i /\ t,r where I is the summation operator. 2. Use 7) with equation 4.33 and n Q at q (the number of reservoirs at the reference flow level) to find the number of reservoirs needed (n t+1 >. and round this value to the nearest integer. 3. If n t+ j ^ n t interpolate the reservoir outflows to compute a new Xjstate vector of dimension n t+ ^. 4. Compute the average flow through each reservoir;
PAGE 170
153 for the first reservoir q = (u +0 J/2 (4.55) t,l k t+1 t,\ and for subsequent reservoirs t,i (,.i + t,ii )/2 <4 56) 5. Compute values for a t an d k t i for the flow levels determined in step 4, using equations 4.31 and 4.35, and the values for a Q and k Q at Q Q 6. Compute parameter a as the average of the delays for the flow through each of the n t+ j reservoirs, as I a J/n i=l C,:L (4.57) 7. Compute the elements of _F^ t a nd Gj( t with equations 4.49 through 4.53. 8. Compute ^ + ^ with equation 1.1. 9. Lag the output by and 10. increment the time and return to step 1. This algorithm allows all parameters of the nonlinear state space model to vary with flow as specified by equations 4.31, 4.33 and 4.35. The thesis of this study is that these relationships can emulate the behavior of the full SaintVenant equations. Results verifying this thesis are presented in Chapter V. The Path to the Final Model The algorithm for the final nonlinear state space routing model performs the necessary sequence of steps to assure that the variations with flow of the moments of the impulse response function for the
PAGE 171
154 linearized SaintVenant equations obey the previously derived theoretical relationships. The approach taken to reach that final algorithm was to begin with a simple algorithm and add complexity only if it improves comparisons with the full unsteady flow model. The major factor complicating the implementation of the nonlinear state space model is the variation with flow of the number of reservoirs, n. There are n elements in the state vector, X_. Letting n vary with flow implies that the dimensionality of X_ changes over time. Initial algorithms were designed to keep n fixed. That simplified the algorithm, but immediately introduced theoretical problems. The lag through a cascade of linear reservoirs is (4.58) Equations 4.31, 4.33 and 4.35, which define the variations with flow of a, n and k, can be substituted into equation 4.58 to produce ( (rl)/m h ~ Q '~ + Â„(m2+r)/ (lm)/m + _0_____ (4 59) Q (lm)/m + Q (lm)/m (4>6Q) Both terms on the right hand side of equation 4.60 vary with flow to the (lm)/m power. If equation 4.33 is simplified as n = n (^61) where n = a constant, then, in order to maintain the variation with flow of the lag indicated by equation 4.60, the variation of k must become
PAGE 172
155 k ~ Q (ml)/m (4.62) Equations 4.62 and 4.35 are not consistent. Holding the number of reservoirs constant warps the variation with flow of k. Because k is also the inverse of K in equation 4.23, the storage computed for each reservoir is no longer that predicted by the theory. If the pure delay, a, is held constant, equation 4.60 becomes t Â£ = Q + Q ^(lm)/m (4.63) which is inconsistent with previous results. The variation of parameter k cannot obey both equations 4.35 and 4.63. In spite of these theoretical shortcomings, several simpler algorithms were developed for a nonlinear state space routing model. In each case, the simplification led to a theoretical inconsistency. The initial algorithms, with n = n Q are: A) 1. Fix a = a at a reference flow level. 2. Compute the lag for each reservoir with equation 4.7. 3. Solve t.i (4.64) (t, 't.i 4. Compute j^ t and _G_x t using the k t j_'s and equations 4.49 through 4.53. 5. Predict _Xt+l with equation 1.1, and 6. Delay the output by a$ Â•
PAGE 173
156 B) 1. Fix k = 1c, at a reference flow level (Compute J_ and G_ for k = k ). 2. Compute the lag for the average flow through the cascade with equation 4.7. 3. Solve n a = t (4.65) t JL k t o 4. Predict X r+1 with equation 1.1, and 5. Delay the output by a t C) 1. Compute k t f r each reservoir with equation 4.7. 2. Compute the lag, to for the average flow through the *t cascade with equation 4.7. 3. Solve n t H 1 r(4 66) C *t i=l t,i 4. Compute _Fy t and ^ t using the k t ^'s and equations 4.49 through 4.53 5. Predict 2^. +1 with equation 1.1, and 6. Delay the output by a t Â• Theoretical inconsistencies caused those preliminary algorithms to emulate poorly the behavior of the primary model. Parameter variations with flow were warped by the improper governing equations. The next step on the path to the final algorithm was to let the number of reservoirs vary with flow. That allowed theoretical consistency to be maintained, but incresed the complexity of the
PAGE 174
157 algorithm. At this point the number of state variables could change over time; therefore, some mechanism had to be developed to produce a new state vector Xjwith n t +i states from an existing Xjwith n t states. All the algorithms studied thus far had assumed that the states (the elements of vector X) were reservoir storages. Reservoir storage and outflow are interchangable, given k t if with equation 4.36. However, when the number of states can vary, storage is an inconvenient choice for the state variable. The approach taken to change the number of states when the states are storage is to maintain continuity. That is, to keep the same total volume of water in the n t+ ^ reservoirs as is originally in the n t reservoirs. Letting n vary when the states are storage created problems because the value of k t> i changes smoothly, while n t+1 is an abrupt change from ru. The simulated outflow for the channel (i.e., from the n t+ reservoir), computed with equation 4.36, has discontinuities when n changes. Those jumps in the flow become more exaggerated as n becomes small, because the percentage change between n t and n t+ ^ increases. To circumvent those problems, the algorithm was rewritten with reservoir outflow as the states. When n changes, the flows are linearly interpolated to produce the new states, and no abrupt changes occur because the downstream outflow is fixed during the interpolation. The algorithm is now essentially in final form. Interpretation of the flow levels at which to compute the parameters, n, a and k, is the last step. Originally, the parameters of the nonlinear state space model varied with the reservoir outflow. The new number of states, n t +l, was found with equation 4.33 using
PAGE 175
v. J, t,i 158 (4.67) Th e a t i 's and k t j's were found with equations 4.31 and 4.35 using =0 (4.68) t,i t,i Equations 4.67 and 4.68 look only at the reservoir outflows to solve for a, n and k. The equations contained in the final algorithm (equations 4.54 through 4.56) use an average of the inflow and outflow to better represent the conditions in the reservoir. These steps lead, finally, to the nonlinear state space routing algorithm presented. All the parameters vary with flow level. In this way the properties of the SaintVenant equations' impulse response function are matched by the nonlinear state space model. Summary A nonlinear state space model that can emulate the behavior of the equations of unsteady flow has been developed. Results obtained by simulating with both the nonlinear state space model and the DWOPER model's numerical solution of the SaintVenant equations are presented in Chapter V. In general, however, any surface runoff or channel routing model can be chosen as the primary model. The strategy employed to emulate any primary model is shown in Figure 4.12. The steps needed to fulfill this strategy are presented in Chapter VI.
PAGE 176
159 SELECT A PRIMARY MODEL i DETERMINE IMPULSE REPONSE PROPERTIES AS FUNCTION OF FLOW LEVEL ASSUME STRUCTURE OF STATE SPACE MODEL i DETERMINE IMPULSE REPONSE PROPERTIES AS FUNCTION OF FLOW LEVEL MATCH IMPULSE RESPONSE FUNCTION PROPERTIES JL NONLINEAR STATE SPACE MODEL AS SURROGATE FOR PRIMARY MODEL Figure 4.12 General Strategy to Emulate Any Primary Model
PAGE 177
CHAPTER V LIMITATIONS OF THE NONLINEAR STATE SPACE APPROACH Introduction The nonlinear state space model developed in Chapter IV models the stationary nonlinear behavior of the full SaintVenant equations with a nonstationary linear structure. Figure 5.1 shows the process used to verify that the nonlinear state space model can emulate the equations of unsteady flow. Given a prismatic channel, the value of m can be found with equation 2.28. For prismatic channel shapes defined by equation 2.23, the value of m will be the same for all flows. However, in general, m can vary with flow if the channel shape if irregular. This possibility is explored in Chapter VI. For results presented in this chapter, m is constant for all flows and can be found directly from equation 2.28. Given m, the variations with flow of parameters a, n and k are specified by equations 4.31, 4.33 and 4.35. Reference points to fix the parameter values for a particular channel shape, flow level and distance from the channel inlet are given in Table 3.4. Recall that these parameter values were computed from the moments of the inflow and outflow hydrographs by first using equations 2.53 and 2.54 to find the mean and variance of the system response function, and then applying equations 2.78 through 2.80. Results from the state space and OWOPER models for large inflow transients are compared in the next section. 160
PAGE 178
161 GIVEN PRISMRTIC CHRNNEL DETERMINE m FOR VRRIOUS FLOH RRNGES BY SIMULATION WITH THE PRIMRRY MODEL i DETERMINE VRRIRTIONS OF a, n RND k WITH FLOW LEVEL INFLOW HYDR0GRRPH5 \ LZ NUMERICRL SOLUTION OF THE SRINTVENRNT EQUATIONS (DWOPER) X FORM NONLINERR STRTE SPRCE MODEL WITH PRRRMETERS a, n AND k Figure 5.1 Verification Process for the Nonlinear State Space Model
PAGE 179
162 Verification of the Nonlinear State Space Model Introduction The nonlinear state space and DWOPER models are compared for various inflow hydrographs in rectangular and triangular channels described in Chapter III. The remaining figures in this chapter are organized as shown in Table 5.1. Each of the 16 figures used to verify the nonlinear state space model corresponds to a figure in Chapter III, comparing the linear state space and DWOPER models. The pairs of corresponding figures are listed in Table 5.2. Results for several singlepeaked inflow hydrographs and a doublepeaked hydrograph are discussed in the following sections. In addition, the issues of continuity, the number of reservoirs required by the state space model, and the CPU times to simulate these results are presented. SinglePeaked Inflow Hydrographs A very rapid transient (
PAGE 180
163 Table 5.1 Organization of Figures 5.2 through 5.17, Comparing the Nonlinear State Space Model with the Full SaintVenant Equations for Large Inflow Transients Rectangular Channel L = 100 miles L = 400 miles o in = 3 hours a. = 6 hours o^ n = 12 hours a. = 24 hours Figure 5.2 Figure 5.4 Figure 5.6 Figure 5.8 Figure 5.3 Figure 5.5 Figure 5.7 Figure 5.9 Triangular Channel in = 6 hours = 12 hours L = 100 miles Figure 5.10 Figure 5.12 L = 400 miles Figure 5.11 Figure 5.13 DoublePeaked Hydrograph = 100000 cfs Q Q = 10000 cfs Q Q = 200000 cfs L = 100 miles Figure 5.14 L = 400 miles Figure 5.15 Figure 5.16 Figure 5.17
PAGE 181
164 Table 5.2 Pairs of Figures with Identical Input Data Comparing the Linear and Nonlinear State Space Models with the Full SaintVenant Equations for Large Inflow Transients Figure Comparing Linear Figure Comparing Nonlinear State Space Model and Full State Space Model and Full SaintVenant Equations SaintVenant Equations 3.28 5.2 3.29 5.3 3.30 5.4 3.31 5.5 3.32 5.6 3.33 5.7 3.34 5.8 3.35 5.9 3.36 5.10 3.37 5.11 3.38 5.12 3.39 5.13 3.40 5.14 3.41 5.15 3.42 5.16 3.43 5.17
PAGE 182
165 O CM cn CD a S ^ OS CD
PAGE 183
166 peak flow a short distance from the channel inlet is significantly lower that the inflow peak of 210000 cfs. This caused problems for the linear state space model when its parameters were selected for Q Q = 100000 cfs, because the flow in the channel was below this level the vast majority of the time. As we have seen in Chapter IV, channel lag decreases as the flow level increases. This means that, since the actual flow was almost always less than the reference flow for the linear state space model, it did not delay the hydrograph sufficiently. Therefore, the entire outflow hydrograph simulated by the linear state space model in Figure 3.28 was earlier than the DWOPER model results. The nonlinear state space model does not adequately lag the rising limb of the inflow hydrograph for a related but different reason. The parameters of the nonlinear state space model are functions of the average flow level in the channel. The average flow in the first reservoir is defined by equation 4.55 as When an inflow hydrograph rises very steeply, this equation gives false information about the actual average flow for that reservoir. The assumption implicit in equation 4.55 is that the average flow is equal to the average of the reservior inflow and outflow. If the inflow disperses very rapidly, equation 4.55 provides a poor representation of the actual average flow. The average flow predicted by equation 4.55 will always be high for the rising limb of a hydrograph. For a very rapidly rising hydrograph it will be significantly high. For a short channel length the problem is compounded because the total number of reservoirs is small and, therefore, the first reservior has a larger
PAGE 184
167 impact on the outflow hydrograph. The parameters of the first reservoir of the nonlinear state space model are computed at a reference flow level higher than actually in the channel. Therefore, as with the linear state space model, the outflow is earlier than the DWOPER model results There is, however, a substantial improvement in the recession, or falling limb, of the state space model outflow hydrograph from Figure 3.28 to Figure 5.2. Once the rapidly varying portion of the inflow hydrograph has entered the channel, the equations for choosing reference flow levels are much more representative of the actual channel conditions. Since the lag and dispersive properties of the nonlinear state space model can vary by changing parameters values, the model is able to match the recession portion of the DWOPER model outflow hydrograph. At a distance of 400 miles from the channel inlet, shown in Figure 5.3, results from the nonlinear state space model are greatly improved over those of the linear state space model in Figure 3.29. The outflow hydrograph from the linear model arrives much too early because its parameters are fixed at = 100000 cfs. The nonlinear state space model is much closer to the lag of the DWOPER model because its parameters vary with flow level in the channel. The rising limb of the nonlinear state space model outflow hydrograph is too early in Figure 5.3 for the same reasons given for Figure 5.2. The situation is ameliorated at 400 miles from the channel inlet, however, because the total number of reservoirs is increased and the first reservoir has a proportionally smaller influence on the model outflow.
PAGE 185
168 a CD i CO OH ts. CC oO O) D r8 8 LR 8 h8 o 31 8 8 _o 8 8 S Â§ K g W 8 K (SJO 0001) SMOUinO 1300U ONU M01JNI ss CO CO o 60 an
PAGE 186
169 Figure 5.4 presents results for a less abrupt inflow hydrograph (a in = 6 hours) at L = 100 miles. Again, there are some discrepancies between the nonlinear state space and DWOPER models for the rising limb of the hydrograph which are caused by an improper representation of the average flow in each reservoir and, therefore, an incorrect choice of state space model parameter values. After the abrupt portion of the inflow hydrograph has entered the channel, the nonlinear state space model agrees very closely with the recession portion of the DWOPER model outflow hydrograph. This is a vast improvement over the linear model in Figure 3.30, which displays the almost symmetrical shape typical of a linear model. A a. = 6 hours inflow hydrograph routed 400 miles in a rectangular channel is shown in Figure 5.5. The timing of the peak predicted by the state space model is very close to that from the DWOPER model. Although the rising limb arrives at the channel outlet too soon, the shape of the recession is a significant improvement over the linear model in Figure 3.31. A comparison of Figures 3.32 and 5.6 provides a classic example of a linear versus a nonlinear model response. The nonlinear state space model for o in = 12 hours at 100 miles in Figure 5.6 produces an outflow hydrograph with sharply rising and slowly falling limbs, and matches the DWOPER model response very well. Outflow from the linear state space model in Figure 3.32 rises too slowly and falls to rapidly. It is doomed by its linear structure to fail when trying to emulate a nonlinear model. At L = 400 miles, in Figures 3.33 and 5.7, the linear versus nonlinear effect is even more pronounced. Although the outflow from the
PAGE 187
170 (sjo ooon SMouino T3Q0W am mum
PAGE 188
171 Q_ az CD Q S O
PAGE 189
172 3: Q_ cc a: CD Se CD Â£ o L
PAGE 190
173 Si o:j ma 9 a:;uu fi QDO x:x:io. tn Cd ZIH Â•Â•alto O i rx Â§ to I Â£*: a. c h ._ I rÂ§ K a . as 8Â§ 3Z 8 W 8 R 8 8 CO N " "Â• Â• Â• in 8 N (sjo ooo n SMOiJino naaow qnu motjni
PAGE 191
174 nonlinear state space model begins to rise too soon, the steepness of its ascent is, in general, equal to that of the DWOPER model. This is in sharp contrast with the Gaussian shape of the linear model in Figure 3.33. The linear state space model performed adequately for the more slowly varying inflow hydrograph (o in = 24 hours) at 100 and 400 miles from the channel inlet as was seen in Figures 3.34 and 3.35. As shown in Figure 5.8, the major improvements introduced by the nonlinear state space model at 100 miles are at the beginning of the rising and end of the recession limbs of the outflow hydrograph. The DWOPER model is matched very well by the nonlinear state space model. At 400 miles from the channel inlet, shown in Figure 5.9, the nonlinear state space model has produced an outflow hydrograph closer to that from the DWOPER model than was possible with the linear model. As mentioned in Chapter III, the incremental improvement from using a linear rather than a nonlinear model would be expected to be greater for a rectangular channel than for a triangular channel. The behavior of a triangular channel is more linear than for a rectangular channel. However, even for a triangular channel, the nonlinear state space model shows improvement over the linear state space model. Comparison of Figures 3.36 and 5.10 for a^ n = f, hours in a triangular channel at L = 100 miles shows that the nonlinear state space model approximates the recession portion of the DWOPER model outflow hydrograph better than the linear model. At 400 miles from the channel inlet, as shown in Figure 5.11, the nonlinear state space model outflow closely resembles the DWOPER model results except for an early rising limb. The timing of the peak flow is much improved over the linear model outflow seen in Figure 3.37.
PAGE 192
175 Q_: en Qii OD: 7 ct: uj 3 S3 :q_ i CO CO ^ I *< = CiJ o I ^ Q C i o I CD CM a v r5 ^ CE en cr a. nS a o a CO o .8 y a a i i CO (SJO 0001) SMOUinO 1300W QNU M01JNI
PAGE 193
176 8 CO o nig o o^ a a Ls s a? o I 9 72 8 (sjo ooon SMOuino i3aow qnu mouni
PAGE 194
177 x! est 9 Sid Â— 3 3 a: [a. u_io IÂ— i.Q t Q O n Â£ a, to u
PAGE 195
178 CD en a> Q S ^ or a >t 3: CD LJ
PAGE 196
179 For 04 = 12 hours in a triangular channel, results from the nonlinear state space model compare well with the DWOPER model, as shown in Figures 5.12 and 5.13. Although the timing of the linear state space model, seen in Figures 3.38 and 3.39, was good, the nonlinearities in the rising and falling limbs are modelled significantly better by the nonlinear state space model. Typical hydrograph durations for natural channels vary widely (Chow 1964, p. 255). An empirical relationshop for the time base of a hydrograph was developed by Snyder (1938) as t T = 3 + (5.1) where t is the lag time to the peak of a unit hydrograph in hours, and T has units of days. The unit hydrograph peak lag time can be expressed as LL t = C (5.2) P C /I where C t s an empirical coefficient (usually equal to 0.35 for a valley drainage area, L is the river mileage through the basin, L is the river mileage from the basin outlet to the basin center of mass, s is the basin slope in feet per mile, and n is an empirical coefficient equal to 0.38 (Linsley, Kohler and Paulhus 1958, p. 207). An approximation of the hydrograph duration can be found with equations 5.1 and 5.2. Hydrographs with o^ n > 12 hours would be expected for most basins. For inflow hydrographs of this duration equation 4.55
PAGE 197
180 Q_ CD Q S CD ^
PAGE 198
181 a. cc en ca 5 s LJ en
PAGE 199
182 adequately captures the rising limb, and the nonlinear state space model emulates well the primary DWOPER model. A DoublePeaked Inflow Hydrograph At 100 miles from the channel inlet, the nonlinear state space model approximates well the DWOPER model for the doublepeaked inflow hydrograph shown in Figure 5.14. As with results presented above, the major improvement of the nonlinear model over the linear state space model is in the shape of the rising limb and in the recession of the outflow hydrograph. In Figure 5.15, a doublepeaked hydrograph is routed 400 miles with the nonlinear state space and DWOPER models. Results from the two models resemble each other very well for the entire outflow hydrograph duration. Figures 5.16 and 5.17 look identical to Figure 5.15. Nonlinear state space model parameters were chosen at Q = 10000 and 200000 cfs in Figures 5.16 and 5.17, respectively. These figures should be compared with Figures 3.42 and 3.42 for the linear state space model. The reference flow level at which parameters are selected is critical to the linear state space model behavior. The variations with flow level expressed in equations 4.31, 4.33 and 4.35 make the nonlinear state space model insensitive to the flow level at which its reference parameter values are chosen. Continuity The linear state space model was guaranteed to preserve continuity (i.e., to make certain that the water leaving the channel equaled the
PAGE 200
183
PAGE 201
184 (SJG 0001) SMOiJinO 13QCW QNU M01JNI
PAGE 202
185 (SJO 0001) SMOIdinO 13Q0W QNU M01JNI
PAGE 203
186 (sjo oooi) SMOuino isaow qnu mouni
PAGE 204
187 water entering). With constant parameters, each reservoir in the linear cascade let all the water that entered pass through with no increase or decrease in quantity. However, there is no such guarantee with the nonlinear state space model. In fact, it is almost certain that the volume of an inflow transient will be modified in some way because of the variations in parameter values. As with any discreteintime model, the nonlinear state space model produces an outflow value each time step. The pure delay, a, changes with the flow level in the channel. The time for which the outflow applies is equal to the time steps taken since the simulation began plus an adjustment for the pure delay. This means that the time period for which an outflow value applies may not be one time step after the previous value simulated. When flow in the channel is increasing, the value of a will decrease. This allows the nonlinear state space model to approximate the steep rising limb of the outflow hydrograph of the DWOPER model. Similarly, the slow decay of the hydrograph' s recession is modelled by increasing the pure delay as flow in the channel decreases. However, this means that the outflow from the last reservoir may apply for less or more than a single time step. During the course of an inflow transient and return to baseflow level, the errors in continuity should compensate. Results presented in Table 5.3 indicate the degree to which this is true. The steeply rising inflow hydrographs (o in = 3 hours) cause the nonlinear state space model to diverge from continuity by 7 to 9 percent. For all other inflows, the nonlinear state space model is within plus or minus 3.3 percent of maintaining continuity.
PAGE 205
Table 5.3 Continuity Check for the Nonlinear State Space and DWOPER Models for Results of Figures 5.2 through 5.17 Figure
PAGE 206
189 The Number of Reservoirs Conceptually, the DWOPER model Is similar to a cascade of nonlinear reservoirs. Each of the distance steps along the channel at which the SaintVenant equations are solved can be thought of as a nonlinear reservoir governed not by equation 4.23, but by the equations of unsteady flow. The number of distance steps used to simulate the results presented above ranged from 10 for a^ n = 24 hours to 160 for a in = ^ hours. The number of reservoirs in the nonlinear state space model is a function of flow level in the channel. It must be an integer each time step, but varies over time as the inflow transient rises and ebbs. The average numbers of reservoirs in the nonlinear state space model cascade are presented in Table 5.4. The ratio of DWOPER model distance steps to nonlinear state space model reservoirs varies from 1:1 to 13:1. We would expect the average number of reservoirs to decrease as the o! value increases because the flow level in the channel is higher for a longer period of time and, therefore, the number of reservoirs decreases. This is not seen in Table 5.4 because the total simulation period was increased as a^ n increased to capture the entire outflow transient. CPU Time Although differences in the nonlinear state space and DWOPER model results shown in Figures 5.2 through 5.17 are generally small, differences in the CPU time required are sometimes large. The CPU times required by all the nonlinear state space and DWOPER model sumulations are presented in Table 5.5. The DWOPER model requires approximately 2
PAGE 207
190 Table 5.4 The Average Number of Reservoirs Used by the Nonlinear State Space Model in Figures 5.2 through 5.17 Rectangular Channel a, = 3 hours in = 6 hours a in = ^ hours a in = 24 hours L = 100 miles 3.533 3.078 3.085 3.201 L = 400 miles 12.212 10.375 11.050 11.894 Triangular Channel L = 100 miles L = 400 miles 'in 2.712 2.413 9.571 8.534 DoublePeaked Hydrograph Q = 100000 cfs Q = 10000 cfs QÂ„ = 200000 cfs L = 100 miles 3.057 L = 400 miles 11.227 11.140 11.243
PAGE 208
191 Table 5.5 CPU Times for the Nonlinear State Space and DWOPER Models on a PRIME 7 50 Computer for Results of Figures 5.2 through 5.17 Rectangular Channel
PAGE 209
192 orders of magnitude more CPU time than does the nonlinear state space model. Summary Results from the nonlinear state space and DWOPER models have been compared for large inflow transients in rectangular and triangular prismatic channels. Inflow hydrographs of varying durations, and a doublepeaked inflow were routed through channels 100 and 400 miles long. The results presented in this chapter have demonstrated that: 1. Many of the important nonlinearities in the full SaintVenant equations are emulated by the nonlinear state space model. 2. The equations for computing average flow in a reservoir generate flow levels that are too high for steeply rising inflow hydrographs. This leads to improper state space model parameter selection and degraded model performance. 3. The nonlinear state space model approximates well the outflow hydrograph recession of the DWOPER model for all inflow hydrographs presented. 4. The problem in Chapter III of an outflow hydrograph that is too symmetrical is solved with the nonlinear state space model 5. The nonlinear state space model typically preserves continuity within plus or minus 3 percent. 6. The DWOPER model requires approximately 2 orders of magnitude more CPU time than the nonlinear state space model. The SaintVenant equations have been our exemplar model to illustrate the new perspective of routing in hydrology. A general strategy to emulate any primary model is the subject of Chapter VI.
PAGE 210
CHAPTER VI HOW TO EMULATE ANY PRIMARY MODEL Introduction Results presented in Chapter V demonstrate that the the nonlinear state space routing model can emulate the full SaintVenant equations. This is only one application of the general approach to nonlinear modelling. The nonlinear state space model's structure should allow it to imitate the behavior of any surface runoff or channel routing model. This chapter describes, in detail, the steps necessary to emulate any primary model. Figure 6.1. summarizes those steps. Select a Primary Model The obvious first step is to choose a primary model. The primary model is chosen based on the data available to describe the physical system and the computational resources committed to simulate the catchment or channel. Calibrate the Primary Model The primary model must be calibrated to the channel or catchment of interest. In general, the parameters of the nonlinear state space model are determined based on simulation with the primary model. The behavior of the nonlinear state space model emulates the primary model, 193
PAGE 211
194 SELECT fl PRIMARY MODEL I CALIBRATE THE PRIMRRY MODEL TO THE CHRNNEL OR CRTCHMENT I USE THE PRIMARY MODEL TO COMPUTE RESPONSE TO INFLOW TRANSIENTS AT SELECTED FLOW LEVELS I COMPUTE MOMENTS OF THE IMPULSE RESPONSE FUNCTION FOR EACH TRANSIENT i ESTIMATE m, a, n AND k AS FUNCTIONS OF FLOW LEVEL CL I K fiT Q m AND a n ARE THE PARAMETERS OF THE NONLINEAR STATE SPACE MODEL Figure 6.1 How to Emulate Any Primary Model
PAGE 212
195 so any errors in calibration of the primary model will affect the performance of the nonlinear state space model. Determine the Impulse Response Properties of the Primary Model Ideally we would determine the weighting function that the primary model uses to transform input signals into output signals. Unfortunately, it is physically impossible to observe the weighting function of a system. As described in Chapter II, for a linear stationary system the weighting function is the timereversed image of the system impulse response. The premise of this work is that the mathematical behavior of a nonlinear stationary system can be represented by a linear model whose impulse response matches that of the nonlinear system at all levels of flow. The impulse response of a nonlinear system changes with flow level; if it did not, the system would be linear. To emulate the changing impulse response of the nonlinear system, the linear model must be nonstationary (i.e., its parameters and, therefore its impulse response, must vary with flow level) Study of the variations in the mean and variance of the impulse response function described in Chapter IV indicates that within flow ranges where the channel shape is regular, the moments of the impulse response function of the SaintVenant equations change as simple algebraic powers of the kinematic wave parameter, m. Therefore, the value of ra can be found from the variations with baseflow level of the moments of the system impulse response. Equations 2.53 and 2.54 express the mean and variance of the system impulse response function. For the mean or variance, the moment of the impulse response function is equal
PAGE 213
196 to the difference of the moment of the system output and the system input. If small inflow transients (i.e., 10 percent of baseflow) are input to the primary model at several baseflow levels, the means and variances of the inflow and outflow can be determined for each baseflow level. Transient inflow hydrographs can be produced with the Fortran program listed in Appendix B. Inputs to this program are: 1. Q the baseflow level, 2. p the incremental transient flow level, 3. Q the flow level about which the equations are linearized, 4. a ^ n the standard deviation of the hydrograph produced, 5. S n BÂ„ and m, which define the channel geometry, and 6. the time step and total period for which the hydrograph will be generated. In general, the means and variances of the impulse response will vary with flow as shown in Figures 4.3 through 4.6. The value of m can be found from the slope of the system impulse response mean (the lag introduced by the channel) versus discharge relationship in log space with equation 4.7. The slope of the lag versus flow level in log space may change with flow level for a general channel shape. Abrupt changes in channel geometry can introduce nonlinearities not accounted for with the algebraic power relationship of impulse response properties with flow. In fact, for a trapezoidal prismatic channel, the value of m must vary with flow. For shallow flows, a trapezoidal channel behaves like a rectangular channel with an m value of 5/3. As the flow deepens, the trapezoidal channel operates more like a triangular channel and has an m value of 4/3.
PAGE 214
197 The value of m can generally be expressed In a piecewise constant versus flow manner as in Figure 6.2. The flow ranges over which m is constant can be determined by the flow levels where the slope is constant in the log(lag) versus log(discharge) relationship. For the particular application explored in Chapter V, the value of m was constant for all flows and could be obtained given the prismatic channel shape. For any prismatic channel shape that obeys equation 2.23, the value of m can be found from Manning's friction formula. In addition, the mean and variance of the system impulse response are expressed as functions of reference flow level by equations 2.56 and 2.58. Estimate a, n and k as Functions of Flow Level The final pieces of information needed for the nonlinear state space model are reference values of the model parameters ( a^ r^ and kg) at any flow level Q i n each of the constant m ranges. The variation of parameters with flow is specified by the value of m. The single reference point at flow level Q fixes the functions for a, n or k. The values of a a, and k. can be found from the impulse response mean and variance for a small transient input on a baseflow, Q within the range of flows for which m is constant. The reference parameter values are computed with equations 2.78 through 2.80. Parameters of the Nonlinear State Space Model The parameters of the nonlinear state space model are now specified. Given m and a f n and k Q at flow level CL for each constant
PAGE 215
198 CD CD O 75 o w w d 0) o Â•2 "2 CD O CM CO 0) Jh he fin (0)901
PAGE 216
199 m flow range, any inflow transient can be modelled with the Fortran program listed in Appendix C. The behavior of the nonlinear state space model mathematically emulates the primary model.
PAGE 217
CHAPTER VII RELATED PREVIOUS WORK Introduction Many other authors have approached the problem of routing surface runoff or channel flow. This chapter describes pertinent previous work in the field of hydrology. Surface runoff and channel routing models related to the nonlinear state space model by either their common cascade of reservoirs structure or approach of modelling nonlinearities by varying parameter values are described in the next section. Observed flow level versus lag time data which corroborates the theoretical relationships developed in Chapter IV are presented. In addition, the nonlinear state space and MuskingumCunge models are compared with the DWOPER model for a 400 mile long rectangular channel. Related Routing Models Early work with models constructed as cascades of reservoirs is attributed to Nash (L957) and to Kalinin and Miljukov (1958). Many authors followed this work in attempts to model a catchment 's instantaneous unit hydrograph (IUH) (Rao, Oelleur and Sarma 1972; Sauer 1973; Pedersen, Peters and Helweg 1980). The IUH is the common term in the field of hydrology for a system's impulse response function. More recent work with cascades of reservoirs has dealt with the development of a state space model of the Nash cascade (SzollosiNagy 1981). The 200
PAGE 218
201 abovementioned works all deal with linear reservoir models that are strictly applicable only over a small flow range around the reference discharge. Amorocho and Brandstetter (1971), and Napiorkowski and Strupczewski (1981) worked with integral forms of the open channel flow equations to parameterize the properties observed in the IUH when routing unsteady flow. The resultant mathematical forms were cumbersome and not amenable to state space formulation. Keefer and McQuivey (1974) took the approach of splitting the flow range into sections and forming a linear model in each range. Calibration of the multiple linearization model was complicated because an impulse response function had to be found for each flow range. They recognized that the IUH varies with flow, but did not study the form of that variation. Other authors developed simple nonlinear models of surface runoff or channel flow (Singh 1964; Dooge 1967; Prasad 1967; Reed, Johnson and Firth 1975; Singh and Buapeng 1981). These studies of the IUH generally expressed the impulse response properties in terms of basin characteristics through complex regression formulas. Reed, Johnson and Firth comment that further study to determine the form of the relationship between physical factors and variable lag is needed. Several authors let model parameters vary with time rather than flow level. Mandeville and O'Donnell (1973) developed time varying expressions for linear reservoir and linear channel equations. The time varying functions for reservoir and channel parameters are fit to observed data. A model in which the parameters change with the season of the year was developed by Chiu and Bittler (1967). They attempted to model the observed nonlinear phenomena by looking for cycles in flow patterns from year to year.
PAGE 219
202 Mein, Laurenson and McMahon (1974) presented a nonlinear storage reservoir model with the proper structure to simulate the nonlinearities in surface runoff or channel flow. The variation with flow of the nonlinear properties is empirically fit to observed data. Their model cannot easily be put in state space form. Schaake (1965, pp. 7692) recognized that the IUH derived from the rising portion of a large surface runoff event was different than the IUH from the recession, and that the impulse response function varied with the magnitude of an event. He concluded that the differences were caused by nonlinear properties of flow routing for large transients. Because of Schaake 's analysis, Valdes, Fiallo and RodriguezIturbe (1979) studied incremental IUHs and concluded that there were not significant differences in the IUHs computed for the rising and falling limbs of hydrographs for small inflow transients. None of the models discussed above has the important combination of properties found in the nonlinear state space model developed in this dissertation. The analysis conducted in the current work has recognized that 1. the variations with flow of model parameters can be described functionally by relating the parameters to the mean and variance of the incremental impulse response for various flow levels, and 2. a state space structure is invaluable for secondorder analysis by means of estimation theory.
PAGE 220
203 Observed Variations with Flow Level of the System Lag Time The variation of lag with flow level has been known in hydrology for some time. Laurenson (1964) empirically derived a relationship between lag and flow level by analyzing hydrographs from the South Creek. Experimental Catchment near Sydney, New South Wales, for 1956 through 1960. Laurenson plotted the lag versus mean flow values on linear, then semilog, and finally loglog axes. He found that the data essentially describe a straight line in logspace. The data he used and the regression equation he fit to the data are shown in Figure 7.1. The slope of the regression line he fit in logspace is 3.70. From equation 4.22 we obtain a value of m = 1.37 for a slope of 3.70. This value of m falls between those of a rectangular (m = 1.67) and a triangular (m = 1.25) channel. Askew (1970) continued on the same path as Laurenson, recognizing the relationship between lag time and flow level. He fit slightly different data than did Laurenson for South Creek, and determined regression equations for lag versus mean flow for four other catchments in Australia. The regression equation exponents and the associated m values found with equation 4.22 are presented in Table 7.1. All the values for m fall close to or within the 1.25 to 1.67 range expected theoretically for triangular and rectangular cross section shapes. Comparison with the MuskingumCunge Model The MuskingumCunge routing model is presented in Appendix A as equations A. 7 through A. 9. This model is developed from a storage routing equation using kinematic wave theory (Cunge 1969). The
PAGE 221
204 K > 3.70 STORAGE DELAY TIME K (HOURS) LAG t t Figure 7.1 Lag Versus Mean Discharge Relationship, South Creek Experimental Catchment, University of New South Wales
PAGE 222
205 Table 7.1 Values of m for Several Catchments
PAGE 223
206 MuskingumCunge model Is the most versatile and physically relavent of the storage routing models (Fread 1982). The nonlinear state space and MuskingumCunge models are compared with the DWOPER model in Figure 7.2. A singlepeaked inflow transient was routed through a 400 mile long rectangular channel with each of the three models. Parameters of the nonlinear state space model are the same as in Chapter V, for 400 mile long rectangular channels. Parameters of the MuskingumCunge model are derived from channel properties as shown in equations A. 8 and A. 9. The MuskingumCunge model matches, in general, the shape of the outflow hydrograph produced by the DWOPER model. However, results with the MuskingumCunge model peak slightly late and high, and remain above the DWOPER model outflow throughout the recession. The nonlinear state space model is early on the rising limb, but the timing of the peak flow and the recession match the DWOPER model excellently. The ratio of CPU time required by the MuskingumCunge model to that of the nonlinear state space model is approximately 25:1. The distance step used in the numerical solution of the MuskingumCunge equations was 5 miles; implying that the MuskingumCunge model used 27 subreaches to simulate the example channel. The average number of reservoirs needed by the nonlinear state space model was 8.057. This ratio of the number of distance steps to reservoirs of 3:1 is consistant with results presented in Chapter V, and is undoubtedly a contributing factor to the additional CPU time needed by the MuskingumCunge model.
PAGE 224
207
PAGE 225
CHAPTER VIII CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK Conclusions The work presented herein consists of a new perspective of routing in hydrology and its application. The new perspective was developed in Chapter I. A particular primary model was then chosen to demonstrate the applicability of the new perspective. The exemplar primary model was the set of equations of unsteady flow in open channels. General conclusions about the new perspective are: 1. The essence of this new perspective of routing is its twostep approach to the problems encountered with existing routing technology. First, the physics of a system are modelled with a primary model; then the mathematics of the primary model are emulated by a state space model. The state space model is a model of the primary model. 2. Stationary nonlinear behavior of a primary model can be emulated by a nonstationary linear model having a state space formulation. 3. The nonlinearities present in the primary model can be captured by matching, at various reference flow levels, the incremental impulse response properties of the primary model with a linear state space model structure. 4. Using a primary model to account for the observed physics, and mathematically emulating the primary model's behavior with a 208
PAGE 226
209 state space model whose parameters vary with flow level is appropriate for modelling surface runoff and channel routing phenomena for downstream flow in prismatic channels. In addition, several conclusions derived from the application of the new perspective to the SaintVenant equations are: 1. The impulse response function of the linearized SaintVenant equations in a prismatic channel has the form of an inverse Gaussian pdf. 2. For small (10 percent of baseflow) inflow transients, the inverse Gaussian pdf models well the behavior of the full equations of unsteady flow. 3. The parameters of the inverse Gaussian pdf can be computed from the moments of the channel impulse response function, which is determined by simulating with the primary model for small inflow transients. 4. For a prismatic channel, the inverse Gaussian pdf parameters can be found directly as functions of the channel geometry. 5. Parameters of the inverse Gaussian pdf are related in a very simple way to parameters of the gamma pdf. 6. A state space model, based on the concept of a cascade of reservoirs, has parameters identical to the 3parameter gamma pdf. 7. A linear state space model (i.e., one whose parameters are constant with flow level) is an adequate surrogate for the full SaintVenant equations for small inflow transients in a prismatic channel, but its behavior deterioriates when large transients are routed.
PAGE 227
210 8. The mean and variance of the impulse response function of the equations of unsteady flow vary as simple algebraic powers of flow level; the expressions for the exponents are functions of the channel geometry. 9. A nonlinear state space model, whose parameters vary to match the changes with flow level of the impulse response function's mean and variance, is an excellent surrogate for the NWS DWOPER model, which is a numerical solution of the full SaintVenant equations. 10. The DWOPER model, which numerically solves the SaintVenant equations in a very efficient manner, requires approximately 2 orders of magnitude more CPU time on a PRIME 750 computer than does the nonlinear state space model. 11. The ability of the state space model to emulate the full SaintVenant equations deteriorates for very abrupt inflow hydrographs, such as from a sudden dam failure, or a small, steepsloped catchment. 12. The current form of the nonlinear state space model is not suitable for modelling channels with lateral inflow or tributaries; however, this limitation should be easily overcome 13. Backwater conditions and upstream wave movement are not accounted for because of the underlying structure of the state space model, which is based on the diffusion equation.
PAGE 228
211 Recommendations for Future Work Some of the inadequacies of current routing technology listed in Chapter I are: 1. it produces no objective estimate of uncertainty, 2. repetitive computations such as in flood frequency analyses can be costly for fully nonlinear models, and 3. error properties of models and data are not expressly considered. The current work directly addresses only item 2 above. The fully nonlinear state space model developed in this work is several orders of magnitude more CPU time efficient than a numerical solution of the SaintVenant equations. Although the state space model has a structure which is amenable to the application of estimation theory, and could produce an objective estimate of uncertainty and consider error properties, that step remains for future work. Several intermediate tasks are required to formulate a filter theory model from the nonlinear state space model. These tasks include: 1. developing a stage versus discharge (rating curve) model which includes a stochastic component to account for shifts in the rating curve, 2. determining a technique to specify values for the filter parameters, and 3. addressing the theoretical problem of filtering when the number of model states varies with time. The nonlinear state space model routes discharge on a catchment or in a channel. Observations of hydrographs are most often made in terms of stage or depth of flow. The relationship between stage and discharge
PAGE 229
212 is not a simple one (Fread 1973). For a given stage, the discharge is higher as the rising limb of a hydrograph passes a given point than during the recession. This is because, during the rising limb, the water level upstream is relatively higher than downstream, so the forces moving the water along the channel are greater. During a recession, with the wave crest now downstream, the water surface slope is less and, therefore, a smaller force is acting on the water. In addition to the variations in rating curves due to water surface slope, channel bed form changes can have dramatic effects, especially in sand bed rivers (Simons, Stevens and Duke 1973). A suitable rating curve model structure compatible with the nonlinear state space routing model and yet accounting for the uncertainties inherent in the stagedischarge relationship must be developed. Kalman filtering is a technique that combines the information in a model and measurements of a process to provide a better estimate of the current conditions of the process than either the model or measurements could give alone (Kalman 1960). A Kalman filter algorithm can be implemented as equation 1.1 and X = F 'X + G U t+1 X,t t X,t t+1 Y = H *X + V (8.1) _t _t _t t where X, F, G and U_ are defined as in equation 1.1, e{u c } = E{V = L Â— t J E ivWi =v 6 i
PAGE 230
213 in which Y^ is the vector of observations at time t, H f is the measurement matrix, Vj. is the vector of measurement errors, q,. is system noise covariance matrix, Rj. is the measurement noise covariance matrix, T indicates a matrix transpose, E is the expectation operator, and 6 is the Dirac delta function which is equal to 1 when i = and otherwise is zero (Jazwinski 1970, p. 269). Equation 1.1 is the structure with which the nonlinear state space model was developed. In a Kalman filtering algorithm it is called the system equation and is responsible for propagating the system dynamics in time. Equation 8.1 is called the measurement equation and is used to incorporate information from measurements of the state vector, X_, into the filter algorithm. The filter parameters that must be specified are the matrices 0_ and R_. Determining values for Q_ and R_will require a major portion of the effort to form an estimation theory model of surface runoff and channel routing with the nonlinear state space model. A theoretical issue that must be addressed before the nonlinear state space model can be used in a Kalman filter algorithm is how to account for a varying number of states. The technique developed in Chapter IV to interpolate the outflows from each reservoir of letting n vary works well when one is simulating with equation 1.1 alone. When, however, the Kalman filter is processing, it propagates not only the
PAGE 231
214 state vector, but also the covariance matrix of the error in X. The system error covariance matrix has a dimension of n by n, with elements that cannot simply be averaged, as is done with the state vector. The first step in attacking this problem should be to explore the literature of other physical science disciplines where estimation theory has been extensively employed. The nonlinear state space model developed in this dissertation is based on a continuous function analysis. The approach taken was to work with the gamma pdf, a continuous function. Differential equations that match the gamma pdf were developed as the nonlinear state space model structure. This restricted us to small time steps in the actual model implementation since we were trying to emulate a continuous function. An alternative approach is to work with a discrete function for the incremental impulse response function (i.e., the pulse response function). The general strategy would be to find a discrete function with the analytical form of the pulse response, and then equate the function and pulse response by matching their moments. Time steps for the model derived from a discrete analysis would be limited only by the frequency properties of the inflow hydrograph, not by numerical aspects of approximating a continuous function. By specifying the changes with flow level of the impulse response function in terms of its mean and variance, we have limited the shape of impulse response functions that can be approximated by the nonlinear state space model to unimodal ones. A different parameterization of the impulse response function could allow multipeaked system responses to be simulated with a nonstationary linear model structure. An analysis in the spectral domain may be the appropriate approach for this problem.
PAGE 232
215 Surface runoff and channel routing comprise only a portion of the hydrologic cycle (Eagleson 1970, pp. 58). Subsurface flow routing is another critical path in the constant motion of water. The highly nonlinear nature of subsurface flow has been modelled with 2and 3dimensional versions of the equations of motion (Eagleson 1970, pp. 271282). Groundwater flow is modelled with a conceptual structure in the NWS Sacramento Soil Moisture Accounting (SMA) model (Burnash, Ferral and McGuire 1973, pp. 1112). Models in state space form have been developed for the NWS SMA model (Kitanidis and Bras 1978, Goldstein and Larimore 1980). These models retain the inherent structure of the original NWS SMA model in a set of simultaneous differential equations used as the system equation, 1.1. It may be possible to extend the new perspective of routing to account for the highly nonlinear structure of groundwater flow. The approach is not clear, however, since the nonlinear properties do not seem to be simple functions of flow level, but instead are related to the interactions of many components. A final area for future work is, of course, to use the nonlinear state space routing model on real catchments and channels. Work, with real data was eschewed in this dissertation because the intent of the new perspective is to use the state space model as a surrogate for a model of the catchment or channel physics. The objective of using an exemplar primary model was to develop the theoretical basis for relating the primary and state space models. The strategy to employ the new perspective has been presented in Chapter VI. This strategy and the Fortran programs listed in Appendices B and C must be tested for actual catchments and channel reaches.
PAGE 233
APPENDIX A A REVIEW OF SOME HYDROLOGIC ROUTING MODELS The complete SaintVenant equations of gradually varied, unsteady flow in an open channel are B .8y + MAv)_ (A#1) 3 1 dx and 1 9v 3y v 3v q r \ tk \ gTFTxgTx f gA v x ; where t = time, x = distance along the channel, y = depth of flow, v = average velocity of the flow, A = cross sectional area of the flow, B = surface width of the flow, g = the acceleration of gravity, q = lateral inflow per unit length, u = the x component of the velocity of the lateral inflow, Sf = the friction slope, and SÂ„ = the channel bottom slope (Henderson 1966, p. 287). The conservation of mass is described by the continuity equation, A.l. The conservation of momentum is expressed by the equation of motion, A. 2. The terms of equation A. 2, from left to right, are measures of the local acceleration, the pressure force, the convective acceleration, the 216
PAGE 234
217 friction force, gravity and the acceleration of the lateral inflow. The quasilinear, hyperbolic system of equations defined by equations A.l and A. 2 has no known analytical solution. Numerical solutions are possible when initial and boundary conditions are specified. All channel routing models obey the continuity equation, A.l, but may be categorized based on which terms of equation A. 2 they consider. Certain routing models do not consider any terms of equation A. 2. They are generally called storage routing models and are based on the continuity equation alone. Equation A.l is usually rewritten as IOJ* (A.3) At where I is the channel inflow, is the channel outflow, and AS is the change in storage in the channel during a At time interval. The storage, inflow and outflow may be related by S = K'[xI + (1X)*0] (A. 4) where K and X are model parameters. Reservoir routing models, such as those developed by Puis (1928) and Goodrich (1931) set X = in equation A. 4. Storage depends only on outflow, and equation A.3 can be expressed in centered finite difference form as I+I 0+0 SS _J L _J 2. = _2 L (A. 5) 2 2 At where the subscripts 1 and 2 indicate the previous and current time steps, respectively. Equation A. 5 can be solved for the unknowns at the current time step as
PAGE 235
218 2S 2S r?+0 =1 +1 + rL ~ (A. 6) At 2 1 2 At 1 The outflow 0, can be found given an S0 relationship from observed data. The wellknown Muskingum model can be found by substituting equation A. 4 into equation A. 5 to yield = C Â• I + C I +C0 (A. 7) 2 2 11 2 1 where C Q = (K> X At/2)/C C : = (KX + At/2)/C 3 C 2 = (K K'X At/2)/C 3 and C 3 = K K'X + At/2 (Chow 1964, p. 2540). In the original Muskingum model as developed by McCarthy (1938), the parmaters K and X were determined from observed inflowoutflow data. Cunge (1969) developed the Muskingum model from kinematic wave theory with a singlevalued stagedischarge relationship. In the MuskingumCunge model the parameters K and X of equation A. 7 are specified as K = Ax/c (A. 8) and X = L[\ /(B c'S Ax)] (A. 9) 2 where c = Â— is the kinematic wave speed at reference discharge dA v Ax = the reach length, S = the channel bottom slope, and B = the channel width at Q Q
PAGE 236
219 Kalinin and Miljukov (1958) developed another model related to the Muskingum model, which is expressed as 2 = 1 + ^1 1^ K 1 + ^2 X l^ K 2 (A 10) where K = 1 e < c At > /Ax 1 K = 1 K Â•Ax/Cc'At), and 2 1 Ax = Q /S Â• (Ah/AQ) in which Ah/AQ is the slope of the stagedischarge relationship. The Muskingum model is equivalent to equation A. 10 If K = Ax/c and X = (Fread 1982). In the Lag and K model (Linsely, Kohler and Paulhus 1949, pp. 230232) the inflow is first lagged, or delayed, and then attenuated with the expression = [i +1 (1 2K/At)]/(l + 2'K/At) (A. 11) 2 1 2 1 This model is linear if the lag and K factors are constant. The lag and K may vary with flow level to give an empirical nonlinear model. Complete dynamic models retain all the terms of equation A. 2. They provide the most accurate models available, hut must still be considered somewhat empirical in nature, since the effects of varying channel geometry are lumped into an empirical roughness coefficient (Weinmann and Laurenson 1979). In most practical applications the value of S n is several orders of magniture larger than the acceleration terms of equation A.l (Henderson 1966, p. 364). A reasonable approach to simplify the full SaintVenant equations is to neglect the acceleration terms of equation A. 2. Various models which are simplified in this way
PAGE 237
220 are called approximate dynamic, diffusion analogy or kinematic wave models (Weinmann and Laurenson 1979). Some of the similarities among the approximate routing models can be seen by rearranging the momentum equation, A. 2 into a looped rating curve formula. This is accomplished by beginning with a flow resistance formula in general form Q = CAR 1 "/"^ (A.12) where C = an empirical resistance coefficient, R = the hydraulic radius, and m = and empirical coefficient. For steady, uniform flow, Q u we can substitute S Q for S f to give Q = CAR m /S (A. 13) By substituting into equation A. 13 for O A* R m from equation A. 12, for S f from equation A. 2, and ignoring the lateral inflow term we can derive a general looped rating curve formula Q = Q { 1 1 3y v 3v 1 3v ; 1/2 S^ST S g 57 S ^ Ht (A. 14) 1 kinematic wave (i.e., steady, uniform flow) diffusion analogy
PAGE 238
221 Â§., (A.i5) and the specified part of equation A. 14. Combining equation A. 15 and the diffusion analogy terms of equation A. 14 produces the equation for a convective diffusion model 80 3Q 8 2 Q / wn *rt + C'r= D Â— + cq (A. 16) at H3x2 where for regular channels c = tt't^ and B dy D = 2*BS The parameters c and D control the convective and attenuation properties, respectively. D is called the diffusion coefficient. Equation A. 16, without the lateral inflow term, approximates equation 2.35 for small Froude numbers. Equation A. 16 was originally proposed by Hayami (1951) with constant parameters c and D. This is essentially the equation solved by the linear state space model developed in this dissertation. By allowing the parameters c and D to vary, equation A. 16 describes a nonlinear model (Price 1973). This is the concept behind the nonlinear state space model. Both Cunge (1969) and Koussis (1976) have demonstrated that with proper choice of parameters, the Muskingum model is a second order approximation of the diffusion equation. For a kinematic wave model the discharge is always a singlevalued function of depth (i.e., the steady, uniform discharge) because from equation A. 14 Q = Q (A. 17)
PAGE 239
222 The kinematic wave model does not account for wave attenuation. The kinematic wave equation can he derived by letting D = in equation A. 16 and rearranging to give where c is now called the kinematic wave speed. If c is constant, equation A. 18 describes a linear model. A nonlinear kinematic wave model arises if c varies. Linearization of equation A. 2 about a reference flow level leads to several other routing models. A form of the wellknown telegraph equation can be obtained by linearizing the friction term, neglecting the convective acceleration, and assuming zero bottom slope and a rectangular channel, as 9 2 v 3 2 v 3v /. Q >. g.y. = __ + g.C J(A. 19) 3X 2 3t 2 W where C Q s a constant related to the linearized friction term (Fread 1982). Lighthill and Whitham (1955) and Harley (1966) developed linear models of the full SaintVenant equations. The complete linear equation obtained by Harley is 3 2 q ^ 9 2 q d 2 y 3 v2).^JL2v 3Â— Q J 9x 2 ^' 3 t 3t 2 (A. 20) u where q is a unit discharge in a unit width channel, SÂ„ is the bottom slope, and
PAGE 240
223 the SaintVenant equations have been linearized about a reference velocity, v = q Q /y (Harley 1966). Equation A. 20 is a particular case of equation 2.33 for a rectangular channel with m = 3/2. Because Harley used a Chezy friction factor in the derivation of equation A. 20, the value of m = 3/2 is appropriate for a rectangular channel. The impulse response function of equation A. 20 is h(x,t) = e P X '6(tx/C ) + nfx/C x/C >e S x r#t l{2'TTm}/m ^12 (A. 21) where C = v + /g*y 1 C = v /gy 2 F = v //gy P = S (2F)/[2y (f 2 +f)] r = S v (2+F 2 )/(2y Â• F 2 J s = S Q /(2.y o ] m =/ (tx/C )(tx/C ) 1 2 l{'} is a first order Bessel function of the first kind, and fi() is the Dirac delta function. Equation A. 20 is related to the diffusion analogy model of Hayami (Chow 1959, pp. 601604). As with any linear model, the response of the complete linearized model is very dependent on the reference flow level (Bravo, Harley, Perkins and Eagleson 1970, pp. 3539).
PAGE 241
APPENDIX B A FORTRAN PROGRAM TO GENERATE HYDROGRAPHS WITH INVERSE GAUSS UN PDF SHAPES IN A PRISMATIC CHANNEL c c C THIS PROGRAM GENERATES AN INVERSE GAUSSIAN SHAPE HYDROGRAPH C FOR RECTANGULAR (EM5/3) OR TRIANGULAR (EM*4/3) CHANNELS FOR C SPECIFIED WAVE DISPERSION PROPERTIES C c C C THIS PROGRAM RUNS ON A PRIME 750 COMPUTER USING FORTRAN 77 C C UNIT 1 IS THE KEYBOARD FOR AN INTERACTIVE TERMINAL C c C C THE VARIABLES ENTERED IN FREE FORMAT ARE: C C QO THE BASE FLOW LEVEL (CFS) C QP THE INCREMENTAL FLOW ADDED TO THE BASE FLOW (CFS) C RQFACT THE REFERENCE FLOW LEVEL, EXPRESSED AS A DECIMAL C FRACTION OF QP C CT THE STANDARD DEVIATION OF THE HYDROGRAPH GENERATED C (HOURS) C BS IF EM=5/3, TOP WIDTH OF THE CHANNEL (FEET) C IF EM=4/3, SIDE SLOPE OF CHANNEL (RISE/RUN) C CM MANNINGS FRICTION COEFFICIENT C M =0, FOR A RECTANGULAR CHANNEL C =1, FOR A TRIANGULAR CHANNEL C SOW CHANNEL BOTTOM SLOPE (FEET/MILE) C DELHR TIME INTERVAL BETWEEN ORDINATES GENERATED (HOURS) C LHOUR LAST HOUR GENERATED C c C C THE INPUT DATA, INTERMEDIATE VALUES, AND THE OUTPUT C TIMES AND DISCHARGES ARE WRITTEN TO THE FILE DESIGNATED C c C CHARACTER*80 FILENAME DIMENSION H(5000) PRINT *, 'ENTER QO, QP, RQFACT, CT, BS.CM.M, SOM.DELKR, LHOUR' 224
PAGE 242
225 ^EAD(1,*)Q0,QP,RQFACT,CT,BS,CM,M,S0M,DELHR,LH0UR IF(QO.LT.0.)GO TO 11 RSOSOM/5280. IF(M.EQ.l)GO TO 5 C c C C RECTANGULAR CHANNEL EM5./3. RQ(QO+RQFACT*QP)/BS RY(RQ*CM/(1.486*SQRT(RSO)))**(l./EM) RVRQ/RY FNRV/SQRT(32.2*RY) RRF1.FN*FN*(1.EM*(2.EM)) RD0.5*RQ/RSO*RRF GO TO 8 C c C C TRIANGULAR CHANNEL C C c 5 EM4./3. Z=BS RQ=QO+RQFACT*QP RY((RQ*CW*2.**(2./3.))/(1.486*SQRT(RSO)*Z))**(l./(2.*EM)) RV=RQ/(Z*RY*RY) FN=RV/SQRT(32.2*RY) RRF=1 .FN*FN*( 1 .EM*( 2 .EM) ) RD=0 5 *RQ / ( RSO* ( 2 *Z*RY ) ) *RRF C c C C THE FOLLOWING CODE IS COMMON TO BOTH CHANNEL SHAPES C C C ~ """"" ~~ 8 CTS=CT*360O. RRC=EM*RV EXO=(((CTS*RRC)**2)*RRC)/(2.*RD) CC2=4.*RD DXO=EXO*RSO/RY RTP((3.*RD)/RRC**2)*(SQRT(l.+(RRC*EXO/(3.*RD))**2)l.) RHP=EXO/SQRT(4.*3.1416*RD*RTP**3)*EXP((EXORRC*RTP)**2/(CC2*RTP)) RQP=QP/RHP CCl=RQP*EXO/SQRT(4.*3.1416*RD) NSTEPS=LHOUR/DELHR IF(NSTEPS.GT.50O0)NSTEPS5OO0 C c
PAGE 243
226 C C GENERATE INVERSE GAUSSIAN SHAPE HYDROGRAPH, H C C c TIMINODELHR*3600. DO 10 Il.NSTEPS TIMEI*TIMINC H(I)QO+(CCl/(TIME**1.5))*EXP((EXORRC*TIME)**2/(CC2*TIME)) 10 CONTINUE C c C STORE RESULTS IN AN OUTPUT FILE C C PRINT *, 'ENTER FILE NAME FOR DISCHARGES AND TIMES' C READ (1,1) FILENAME 1 FORMAT (A80) OPEN (FILEFILENAME UNIT1 05 ) WRITE(105,500)QO,QP,RQFACT,CT,BS,CM,EM,SOM,DELHR,LHOUR, 1 RQ,SOM,RSO,RY,RV,FN,RRF,RD,CTS,RRC, 2 EXO,CC2,DXO,RTP,RHP,RQP,CCl C 500 FORMAT ('QO', Gil. 4,', QP'G11.4,', RQFACT' ,G11 .4/ 1 'CT',G11.4,', BS',G11.4/ 1 'CM=',G11.4,', EM', Gil. 4,', S0M',G11.4/ 1 'DELHR',G11.4,', LHOUR'.Ill/ 1 'RQ',G11.4,', S0M',G11.4/ 1 'RS0=',G11.4,', RY',G11.4/ 1 'RV=',G11.4,', FN', Gil. 4/ 1 'RRF',G11.4,', RD',G11.4/ 1 'CTS=',G11.4,', RRC=',G11.4/ 1 'EX0=',G11.4,', CC2',G11.4/ 1 'DX0=',G11.4,', RTP',G11.4/ 1 'RHP', Gil .4,', RQP=',G11.4/ 1 'CC1=',G11.4) C WRITE(105,501)(H(I),I1,NSTEPS) 501 FORKAT(8F10.2) C 11 STOP END
PAGE 244
APPENDIX C A FORTRAN PROGRAM FOR THE NONLINEAR STATE SPACE ROUTING MODEL C c C THIS PROGRAM SIMULATES WITH THE NONLINEAR STATE SPACE MODEL C c c C THIS PROGRAM RUNS ON A PRIME 750 COMPUTER UNDER FORTRAN 77 C C UNIT 1 IS THE KEYBOARD FOR AN INTERACTIVE TERMINAL C c c C THE INPUT READ FROM A FILE IN FREE FORMAT IS: C C AO THE PURE DELAY AT THE REFERENCE FLOW LEVEL C NO THE NUMBER OF RESERVOIRS AT THE REFERENCE FLOW LEVEL C (THIS MAY BE A NONINTEGER VALUE) C KO THE PARAMETER K AT THE REFERENCE FLOW LEVEL C QO THE REFERENCE FLOW LEVEL C M THE KINEMATIC WAVE PARAMETER C DT THE TIME STEP C NSTI THE NUMBER OF VALUES OF INFLOW ENTERED C (DT*NSTI IS THE PERIOD SIMULATED) C RHO THE TIME WEIGHTING FACTOR (RANGE FROM TO 1) C DEBUG =0, NO DEBUG PRINTOUT, C =1, PRINTOUT FOR EACH RESERVOIR FOR EACH TIME STEP C INFLOW THE INFLOW TIME SERIES (NSTI VALUES) C c C C THE OUTFLOW FROM THE DOWNSTREAM RESERVOIR AND THE CORRESPONDING C TIMES ARE STORED IN ARRAYS RESOUT AND TIMOUT, RESPECTIVELY. CHARACTER*80 FILENAME LOGICAL LINEAR INTEGER*4 TOTRES, DEBUG REAL*4 K M KO NO .KOLOG LGOLOG INFLOW INFBAR K I INFLOG NO LOG C0MM0H/X1X/ F(100,100),G(100),X(100),XP(100),GU(100),K(100), 1 TLAG(100),AX(100),QLIN(100) COMMON/X2X/INFLOW(5001 ) ,RESOUT(5001 ) ,TIMOUT(5001 ) NRD=100 227
PAGE 245
228 HWR101 C PRINT*, 'ENTER INPUT FILE NAME' C READ (1,1) FILENAME 1 FORMAT (A80) OPEN( ONITNRD FILEFILENAME ) C PRINT*, 'ENTER OUTPUT FILE NAME' C READ(1,1)FILENAME OPEN(UNITNWR, FILEFILENAME) C READ(NRD,*)AO, NO, K0,Q0,M,DT,NSTI,RHO, DEBUG IF(N0.EQ.0.)GO TO 999 C READ(NRD,*)(INFL0W(I),I1,NSTI) C NMAX100 LINEAR. FALSE. IF(M.EQ.1.)LINEAR.TRUE. IF(LINEAR)GO TO 400 TOTLGOAO+NO/KO BLAG(1.M)/M C C P = 2/3 FOR MANNINGS FRICTION FORMULA C R IS A FUNCTION OF M FROM THE PRISMATIC SHAPE FORMULA C P2./3. R(l.+PM)/P BN=(R1.)/M BK(M2.+R)/M QOLOGLOGIO(QO) KOLOG=LOG10(KO) AOLOG=LOG10(AO) N0LOG=LOG10(N0) LGOLOG=LOG10(TOTLGO) C INFLOG=LOG10 ( INFLOW( 1 ) ) TOTLAG=10.**(BLAG*(INFLOGQOLOG)+LGOLOG) N=10.**(BN*(INFLOGQ0LOG)+N0LOG) + 0.5 IF ( N GT NMAX ) NNMAX IF(N.LT.1)N=1 KI=10.**(BK*(INFLOGQ0LOG)+K0LOG) AI=TOTLAGNI/KI GO TO 401 C 400 KI=K0 N=N0+0.5 A=A0 C 401 S0=INFLOW(l) INFBAR=INFL0W(1) V0LIN=0.
PAGE 246
229 VOLOUTLR0. TOTRES0 C DO 2 11,11 K(I)KI 2 X(I)S0 C DO 10 Jl.NSTI C IF(J.GT.l)INFBAR(INFLOW(J)+INFLOW(Jl))/2. C IF (LINEAR) GO TO 405 C C COMPUTE AVERAGE OUTFLOW FROM ALL RESERVOIRS C AS FLOW TO DETERMINE NEW NUMBER OF RESERVOIRS C QOBARINFBAR DO 705 11,1 705 QOBAR"QOBAR+X(I) QOBARQOBAR/(N+l) C NOLDN C C COMPUTE NEW NUMBER OF RESERVOIRS BASED ON QOBAR C N10.**(BN*(LOG10(QOBAR)QOLOG)+NOLOG) + 0.5 IF(N.GT.NMAX)NNMAX IF(N.LT.1)N1 C TOTRES=TOTRES+N IF(N.EQ.NOLD)GO TO 730 C C IF NUMBER OF RESERVOIRS HAS CHANGED C INTERPOLATE OLD STATES TO GET DISCHARGES FOR NEW C NUMBER OF RESERVOIRS C JJ=0 LCM=N*NOLD C DO 710 I1,LCM IF((l/NOLD)*NOLD.LT.I)GO TO 710 I1=(I1)/N + 1 I2=(I1)/N0LD + 1 QLEFT=INFBAR IF(I1 .GT. 1 )QLEFT=X(I11 ) QRIGHT=X(I1) QDIFF=QRIGHTQLEFT XLEFT=N*(I11) XRIGHT=N*I1 XWANT=NOLD*I2 JJ=JJ+1 QLIN(JJ)=QLEFT+QDIFF*(XWANTXLEFT)/(XRIGHTXLEFT) 710 CONTINUE
PAGE 247
230 C FILL X ARRAY WITH OUTFLOWS FROM EACH NEW RESERVOIR C DO 720 Il.N 720 X(I)QLINU) C 730 ATOT0. DO 3 11,1 C IF(I.EQ.l)QBAR(lNFBAR+X(l))/2. IF(I.GT.l)QBAR(X(Il)+X(I))/2. XXXLOG10 (QBAR)QOLOG AX ( I ) 1 ** ( BLAG*XXX+A0LOG ) K ( I ) 1 .** ( BK*XXX+K0L0G ) TLAG ( I ) 1 ** ( BLAG*XXX+LG0LOG ) ATOTATOT+AX(I) 3 CONTINUE AATOT/N C 405 IF(.NOT.LINEAR.OR.J.EQ.l)CALL FILLFG(F,G,NMAX,N,K,DT,RHO,NWR) C CALL LDMULT(F,X,XP,N,1,NMAX,1) DO 5 Il.N 5 GU(I)G(I)*INFBAR CALL ADD(XP,GU,X,N,1,NMAX,1) RESOUT(J)X(N) NDELAY=A/DT+0.5 ITMOUT=J+NDELAY TIMOUT(J)=ITMOUT*DT IF(DEBUG.EQ.0)GO TO 410 WRITE(NWR,653)J,A,Q0BAR,N,INFBAR,RES0UT(J),TIM0UT(J) 653 FORMATC ><><><><><>< TIME STEP', 15,', A',F10.3,', QOBAR', 1 F10.4,', N=',I4,', INFBAR',F10.3,', RESOUT(J)', 2 F10.4,', TIMOUT(J)',F8.2) WRITE(NWR,654)(X(I) ,11 ,N) 654 FORMAT (' DISCHARG' ,10(1X,F10 .4)/( 10X,10(1X,F10.4) ) ) IF(LINEAR)GO TO 410 WRITE(KWR,655)(K(I),I1,N) 655 FORMATC K',10(1X,F10.4)/(10X,10(1X,F10.4))) WRITE(NVR,648)(TLAG(I),I=1,N) 648 FORMATC TLAG=' ,10( 1X,F10.4)/(10X, 10( 1X.F10.4) )) WRITE(KWR,673)(AX(I),I=1 ,N) 673 FORMATC AX=' ,10( 1X.F10 .4)/( 10X, 10( 1X.F10.4) ) ) IF(N.NE.N0LD)WR1TE(NWR,674)(QLIN(I),I1,N) 674 FORMATC QLIN' ,10( lX,F10.4)/( 10X, 10( 1X.F10.4))) C 410 VOLIN=VOLIN+INFLOW(J)*DT VOLOUTLR=VOL0UTLR+RESOUT(J)*DT
PAGE 248
231 C 10 CONTINUE C AVGRESFLOAT(TOTRES)/NSTI WRITE(NWR,657)VOLIN,VOLOUTLR,AVGRES 657 FORMATC THE TOTAL INFLOW '.G17.6/ 2 THE TOTAL OUTFLOW FROM RESERVOIRS '.G17.6/ 3 THE AVERAGE NUMBER OF RESERVOIRS '.F8.3) C 999 STOP END c c SUBROUTINE FILLFG(F,G,MAX,N,K,DT,RHO,NWR) C C c C FILL F AND G MATRICES FOR X(T+1) F*X(T) + G*U EQUATION C FOR STATES OF DISCHARGE C REAL*4 K DIMENSION F(MAX,MAX),G(MAX),K(MAX) C C FILL MATRIX F(l,J) IROWS, JCOLS C GAMMA1 .RHO C C FILL FROM 1ST TO NTH ROWS C SINCE K IS CONSTANT FOR A ROW C DO 10 Il.N SQUIGL=1 ./K(I)+DT*RHO IF(I.EQ.N)GO TO 6 L=I+1 C C ZERO TO THE RIGHT OF MAIN DIAGONAL C DO 5 J=L,N 5 F(I,J)=0. C C FILL IN VALUE ON MAIN DIAGONAL C 6 F(I,I)=(1./K(I)DT*GAMMA)/SQUIGL IF(I.EQ.l)GO TO 9 C C FILL IN LEFT OFF DIAGONAL C L=I1 F(I,L)=(F(L,L)*RHO+GAMMA)*DT/SQUIGL IF(L.EQ.l)GO TO 9 C C FILL TO THE LEFT OF LEFT OFF DIAGONAL
PAGE 249
232 C FILL FROM RIGHT TO LEFT C DO 8 JJ2.L JIJJ F(l,J)F(L,J)*DT*RHO/SQUIGL 8 CONTINUE C C FILL IN G VECTOR VALUE C 9 IF(I.EQ.1)G(I)DT/SQUIGL IF(I.GT.1)G(I)G(L)*DT*RH0/SQUIGL 10 CONTINUE C RETURN END SUBROUTINE LDMULT(A,B,C,L,M,LDIM,MDIM) C C MULTIPLY TWO LOWER DIAGIONAL MATRICES C DIMENSION A(LDIM,LDIM),B(LDIM,MDIM),C(LDIM,MDIM) DO 20 Il.L DO 20 Jl.M TEMP0 DO 10 K1,1 10 TEMP=TEMP+A(I,K)*B(K,J) 20 C(I,J)=TEMP RETURN END c C SUBROUTINE ADD(A, B,C,N,M,NDIM,MDIM) C C ADD TWO MATRICES C c C DIMENSION A(NDIM.MDIM) ,B(NDIM,MDIM) ,C(NDIM,MDIM) DO 10 1=1, N DO 10 J=1,M 10 C(I,J)=A(I,J)+B(I,J) C RETURN END
PAGE 250
BIBLIOGRAPHY Abramowitz, M., and I. A. Stegun, editors, Handbook of Mathematical Funct ions with Formulas, Graphs, and Mathematical Tables U.S. Department of Commerce, National Bureau of Standards, Applied Mathematics Series 55, U.S. Government Printing Office, Washington, D.C., June, 1964. Amorocho, J., and A. Brandstetter, "Determination of Nonlinear Functional Response Functions in RainfallRunoff Processes," Water Resources Research Vol. 7, No. 5, pp. 10871101, October, 1971. Askew, A. J., "Derivation of Formulae for Variable Lag Time," Journal of Hydrology, Vol. 10, pp. 225242, 1970. Bansal, M. K. "Dispersion in Natural Streams," Journal of the Hydraulics Division, Vol. 97, No. HY11, pp. 18671886, November, 1971. Bravo, C. A. B., B. M. Harley, F. E. Perkins and P. S. Eagleson, _A Linear Distributed Model of Catchment Runoff, Technical Report No. 12 3, Ralph M. Parsons Hydrodynamics Laboratory, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, June, 1970. Burnash, R. J. C, R. L. Ferral and R. A. McGuire, A Generalized Streamflow Simulation System Conceptual Modeling for Digital Computers U.S. Department of Commerce, National Weather Service, Washington, D.C., and State of California, Department of Water Resources, Sacramento, California, March, 1973. Carslaw, H. S., and J. C. Jaeger, Conduction of Heat in Solids Oxford University Press, London, 1959. Chiu, C.L., and R. P. Bittler, "Linear TimeVarying Model of RainfallRunoff Relation," Water Resources Research, Vol. 5, No. 2, pp. 426437, April, 1969. Chow, V. T., Open Channel Hydraulics, McGrawHill Book Company, New York, 1959. Chow, V. T., editorinchief, Handbook of Applied Hydrology: A Compendium of WaterResources Technology McGrawHill Book Company, New York, 1964. Crank, J., The Mathematics of Diffusion Oxford University Press, London, 1956. Cunge, J. A., "On the Subject of a Flood Propagation Computation Method (Muskingum Method)," Journal of Hydraulics Research, Vol. 7, No. 2, pp. 205230, 1969. 233
PAGE 251
234 Daily, J. W., and D. R. F. Harleman, Fluid Dynamics, AddisonWesley Publishing Company, Inc., Reading, Massachusetts, 1966. Dooge, J. C. I., "A New Approach to Nonlinear Problems in Surface Water Hydrology: Hydrologic Systems with Uniform Nonlinearity," International Association of Scientific Hydrology Publication 76, pp. 409413, 1967. Dooge, J. C. I., Linear Theory of Hydrologic Systems Technical Bulletin No. 1468, Agricultural Research Service, U.S. Department of Agriculture, Washington, D.C., October, 1973. Eagleson, P. S., Dynamic Hydrology McGrawHill Book Company, New York, 1970. Fischer, H. B., "The Mechanics of Dispersion in Natural Streams," Journal of the Hydraulics Division, Vol. 94, No. HY6, pp. 187216, November, 1967. Fread, D. L., A Dynamic Model of StageDischarge Relations Affected by Changing Discharge NOAA Technical Memorandum NWS HYDRO16, U.S. Department of Commerce, National Weather Service, Silver Spring, Maryland, November, 1973. Fread, D. L., Numerical Properties of Implicit FourPoint Finite Difference Equations of Unsteady Flow NOAA Technical Memorandum NWS HYDRO18, U.S. Department of Commerce, National Weather Service, Silver Spring, Maryland, March, 1974. Fread, D. L., "NWS Operational Dynamic Wave Model," in Proceedings of Verification of Mathematical and Physical Models 26th Annual Hydraulics Division Specialty Conference, American Society of Civil Engineers, College Park, Maryland, pp. 455464, 1978. Fread, D. L., "Flood Routing: A Synopsis of Past, Present, and Future Capability," in RainfallRunoff Relationship a part of the Preceedings of the International Symposium on RainfallRunoff Modeling, Mississippi State University, May 1821, 1981, Water Resources Publications, Littleton, Colorado, 1982. Fread, D. L., and G. F. Smith, "Calibration Technique for 1D Unsteady Flow Models," Journal of the Hydraulics Division, Vol. 104, No. HY7, pp. 10271044, July, 1978. Goldstein, J. D., and W. E. Larimore, Applications of Kalman Filtering and Maximum Likelihood Parameter Identif icatior! to Hydrologic Forecasting Technical Report TR14801, prepared by The Analytic Sciences Corporation, Reading, Massachusetts, under Contract No. NA79SAC00668 for the U.S. Department of Commerce, National Oceanic and Atmospheric Administration, National Weather Service, April, 1980.
PAGE 252
235 Goodrich, R. D., "Rapid Calculation of Reservoir Discharge," Civil Engineering Vol. 1, pp. 417418, 1931. Harley, B. M., Linear Response Functions of Uniform Open Channels Thesis, Department of Civil Engineering, University College, Cork, Ireland, April, 1966. Hayami, S., On the Propagation of Flood Waves Bulletin No. 1, Disaster Prevention Research Institute, Kyoto University, Japan, December, 1951. Henderson, F. M., Open Channel Flow The Macmillan Company, New York, 1966. Jazwinski, A. H., Stochastic Processes and Filtering Theory Academic Press, New York, 1970. Johnson, J. R., and D. E. Johnson, Linear Systems Analysis The Ronald Press Company, New York, 1975. Johnson, N. L., and S. Kotz, Continuous Univariate Distributions 1 John Wiley and Sons, New York, 1970. Kalinin, G. P., and P. I. Miljukov, "On the Computation of Unsteady Flow Along the Channels by the Use of Reachtravel Curves," Heterology and Hydrology Vol. 7, pp. 1825, 1958. Kalman, R. E., "A New Approach to Linear Filtering and Prediction Problems," Journal of Basic Engineering Vol. 82D, pp. 3546, March, 1960. Keefer, T. N., and R. S. McOuivey, "Multiple Linearization Flow Routing Model," Journal of the Hydraulics Division, Vol. 100, No. HY7, pp. 10311046, July, 1974. Kitanidis, P. K. and R. L. Bras, Real Time Forecasting of River Flows Technical Report No. 235, Ralph M. Parsons Hydrodynamics Laboratory, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, June, 1968. Koussis, A. D., "An Approximative Dynamic Flood Routing Method," Paper LI in International Symposium on Unsteady Flow in Open Channels University of NewcastleuponTyne, England, April 1215, 1976. Koussis, A. D., "Theoretical Estimations of Flood Routing Parameters," Journal of the Hydraulics Division Vol. 104, No. HY1, pp. 109115, January, 1978. Laurenson, E. M., "A Catchment Storage Model for Runoff Routing," Journal of Hydrology, Vol. 2, pp. 141163, 1964.
PAGE 253
236 Liggett, J. A., and J. A. Cunge, "Numerical Methods of Solution of the Unsteady Flow Equations," Chapter 4 in Unsteady Flow in Open Channels, Volume 1, edited by K. Mahmood and V. Yevjevich, Water Resources Publications, Fort Collins, Colorado, 1975. Lighthill, M. J., and G. B. Whitham, "On Kinematic Waves I. Flood Movement in Long Rivers," Proceedings, Royal Society Vol. 229A, pp. 281316, May, 1955. Linsley, R. K. M. A. Kohler and J. L. H. Paulhus, Hydrology for Engineers McGrawHill Book Company, New York, 1958. Liu, C. L., and J. W. S. Liu, Linear Systems Analysis McGrawHill Book Company, New York, 1975. Mandeville, A. N., and T. O'Donnell, "Introduction of Time Variance to Linear Conceptual Catchment Models," Water Resources Research, Vol. 9, No. 2, pp. 298310, April, 1973. McCarthy, G. T., "The Unit Hydrograph and Flood Routing," presented at Conference of the North Atlantic Division of the U.S. Corps of Engineers, New London, Connecticut, 1938. Mein, R. G., E. M. Laurenson and T. A. McMahon, "Simple Nonlinear Model for Flood Estimation," Journal of the Hydraulics Division, Vol. 100, No. HY11, pp. 15071518, November, 1974. Miller, W. A., and J. A. Cunge, "Simplified Equations of Unsteady Flow," Chapter 5 in Unsteady Flow in Open Channels, Volume 1 edited by K. Mahmood and V. Yevjevich, Water Resources Publications, Fort Collins, Colorado, 1975. Napiorkowski, J. J., and W. G. Strupczewski "The Properties of the Kernels of the Volterra Series Describing Flow Deviations from a Steady State in an Open Channel," Journal of Hydrology Vol. 52, pp. 185198, 1981. Nash, J.E., "The Form of the Instantaneous Unit Hydrograph," Proceedings of the General Assembly, International Association of Scientific Hydrology Vol. 3, pp. 121144, 1957. Nash, J. E., "Systematic Determination of Unit Hydrograph Parameters," Journal of Geophysical Research, Vol. 64, No. 1, pp. 111115, 1959. Ogata, K., State Space Analysis of Control Systems PrenticeHall, Inc., Englewood Cliffs, New Jersey, 1967. Pedersen, J. T., J. C. Peters and 0. J. Helweg, "Hydrographs by Single Linear Reservoir Model," Journal of the Hydraulics Division, Vol. 106, No. HY5, pp. 837852, May, 1980.
PAGE 254
237 Prasad, R., "A Nonlinear Hydrologic System Response Model," Journal of the Hydraulics Division Vol. 93, No. HY4, pp. 201221, July, 1967. Puis, L. G., "Construction of Flood Routing Curves," U.S. 70th Congress, 1st Session, House Document 185, pp. 433442, 1928. Rao, R. A., J. W. Delleur and R. S. P. Sarma, "Conceptual Hydrologic Models for Urbanizing Basins," Journal of the Hydraulics Division, Vol. 98, No. HY7, pp. 12051220, July, 1972. Reed, D. W., P. Johnson and J. M. Firth, "A NonLinear RainfallRunoff Model, Providing for Variable Time Lag," Journal of Hydrology, Vol. 25, pp. 295305, 1975. Sauer, V. B., "UnitResponse Method of OpenChannel Flow Routing," Journal of the Hydraulics Division Vol. 99, No. HY1, pp. 179193, January, 1973. Schaake, J. C, Jr., Synthesis of the Inlet Hydrograph Dissertation, Johns Hopkins University, Baltimore, Maryland, 1965. Schaake, J. C, Jr., "Analytical Solution of Flood Wave Attenuation in Prismatic Channels," in Hydrology, Geomorphology and Climate, J. Valdes, coordinator, Simon Bolivar University, Caracas, pp. 147179, 1980. Simons, D. B., M. A. Stevens and J. H. Duke, "Predicting Stages on SandBed Rivers," Journal of the Waterways, Harbors and Coastal Engineering Division Vol. 99, No. WW2, pp. 231243, May, 1973. Singh, K. P., "Nonlinear Instantaneous UnitHydrograph Theory," Journal of the Hydraulics Division Vol. 90, No. HY2, pp. 313347, March, 1964. Singh, V. P., and S. Buapeng, "A Nonlinear Hydrologic Cascade," Journal of Hydrology, Vol. 51, pp. 283293, 1981. Slade, D. H., editor, Meterology and Atomic Energy, 1968 U.S. Atomic Energy Commission, Division of Technical Information, Publication No. TID24190, July, 1968. Snyder, F. F., "Synthetic Unit Hydrographs ," Transactions American Geophysical Union Vol. 19, Part 1, pp. 447454, 1938. SzollosiNagy, A., "State Space Models of the NashCascade, Kinematic and Diffusion Waves," Water Resources Engineering Lulea, University of Lulea, Sweden, Series A, No. 68, April, 1981. Thomas, J. B., An Introduction to Statistical Communication Theory John Wiley and Sons, Inc., New York, 1969.
PAGE 255
238 Valdes, J. B., Y. Fiallo and I. RodriguezIturbe, "A RainfallRunoff Analysis of the Geomorphologic IUH," Water Resources Research Vol. 15, No. 6, pp. 14211434, December, 1979. Weinraann, P. E., and E. M. Laurenson, "Approximate Flood Routing Methods: A Review," Journal of the Hydraulics Division, Vol. 105, No. HY12, pp. 15211536, December, 1979. Woolhiser, D. A., and J. A. Liggett, "Unsteady, OneDimensional Flow over a Plane the Rising Hydrograph," Water Resources Research Vol. 3, No. 3, pp. 753771, Third Quarter, 1967. Zoch, R. T., "On the Relation Between Rainfall and Runoff, Part I," Monthly Weather Review, Vol. 62, pp. 315322, 1934.
PAGE 256
BIOGRAPHICAL SKETCH George Francis Smith, born on July 13, 1951, in Englewood, New Jersey, attended secondary school in Rergenfield, New Jersey. He graduated from the Massachusetts Institute of Technology with a Bachelor of Science degree in civil engineering in 1973. After completing a Master of Engineering in the environmental engineering program at the University of Florida in 1975, he began work with the National Weather Service Hydrologic Research Laboratory in Silver Spring, Maryland. He and Margaret James Smith have lived in Silver Spring since they were married in 1976, except for the 19781979 academic year when they returned to the University of Florida on a fulltime university assignment from the National Weather Service. 239
PAGE 257
I certify that I have read this studv and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. i / C tiuL Wayne C; Huber, Chairman Professor of Environmental Sciences engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fullv adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. John C. Schaake Adjunct Professor of Environmental Engineering Sciences 1 certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. / J^mes P. Heaney / ''Professor of Environmental Engineering Sciences I certify that I have read this study and that in my ooinion it conforms to acceptable standards of scholarly nresentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Thomas E. Bullock Professor of Electrical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarlv presentation and is fully adequate, in scope and qualitv, as a dissertation for the degree of Doctor of Philosophy. ^ y n Antal Maythay Associate Professor of Management and Administrative Sciences
PAGE 258
This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate Council, and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1983 I "t^JAtJf U" O^t^ Dean, College of Engineering lLhJ QL Dean for Graduate Studies and Research
PAGE 259
3 1262 08666 315 9
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID ESJO4M1YV_EM4LKC INGEST_TIME 20171109T20:45:44Z PACKAGE UF00101079_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E9S9LK3QV_Q1OCVJ INGEST_TIME 20170711T22:17:51Z PACKAGE UF00101079_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES

