Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Optimal allocation of ’BOD’ loadings in a tidal river Jamal, Iqbal Badrudin 1986

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1986_A4_6 J35.pdf [ 6.68MB ]
Metadata
JSON: 831-1.0096653.json
JSON-LD: 831-1.0096653-ld.json
RDF/XML (Pretty): 831-1.0096653-rdf.xml
RDF/JSON: 831-1.0096653-rdf.json
Turtle: 831-1.0096653-turtle.txt
N-Triples: 831-1.0096653-rdf-ntriples.txt
Original Record: 831-1.0096653-source.json
Full Text
831-1.0096653-fulltext.txt
Citation
831-1.0096653.ris

Full Text

OPTIMAL ALLOCATION OF *BOD' LOADINGS IN A TIDAL RIVER By IQBAL BADRUDIN JAMAL B.A.Sc, The University of B r i t i s h Columbia, 1978 M.A.Sc, The University of B r i t i s h Columbia, 1980 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Faculty of Commerce and Business Administration) We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1986 @ Iqbal Badrudin Jamal, 19^ 36 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. I t i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of CpMMSR£E AND BpSl^fe£ Am\MSrQXnc>h\ The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date ABSTRACT A methodology Is presented i n this thesis which addresses the water quality manager's problem of maximizing the biochemical oxygen demand loadings in reaches of a t i d a l river, subject to dissolved oxygen concentration regulation at various compliance levels. The non-linear t i d a l dynamics and BOD-DO processes were ex p l i c i t l y accounted for by numerical, finite-difference models incorporated as equality and inequality constraints of a non-linear programming problem solved by a direct search algorithm. This methodology was applied to the Nicomekl River in Surrey B.C., to investigate policy implications and i t s applicability, using a microcomputer, to the resolution of actual pollution management problems. - i i i -T A B L E OF CONTENTS A b s t r a c t i i T a b l e o f C o n t e n t s i i i L i s t o f T a b l e s v i i L i s t o f F i g u r e s v i i i A c k n o w l e d g e m e n t s x 1 . INTRODUCTION 1 2. LITERATURE SURVEY 4 2 .1 I n t r o d u c t i o n 4 2 .2 W a t e r Q u a l i t y M o d e l s 5 2 . 2 . 1 S t e a d y - S t a t e M o d e l s 6 2 . 2 . 2 Dynamic M o d e l s 12 2 . 2 . 3 E s t i m a t i o n T e c h n i q u e s 18 2 . 3 Sys tems O p t i m i z a t i o n T e c h n i q u e s 20 2 . 3 . 1 L i n e a r Programming M o d e l s 21 2 . 3 . 2 N o n l i n e a r Programming M o d e l s 25 2 . 3 . 3 Dynamic Programming M o d e l s 27 2 . 3 . 4 S t o c h a s t i c Programming M o d e l s 32 2 . 3 . 5 C o n t r o l M e t h o d s 35 2 . 4 Summary 38 - i v -3 . MODEL DEVELOPMENT 39 3 . 1 P r o b l e m F o r m u l a t i o n 39 3 . 1 . 1 O b j e c t i v e F u n c t i o n 39 3 . 1 . 2 C o n s t r a i n t s 40 3 . 2 Water Q u a l i t y M o d e l 44 3 . 2 . 1 H y d r o d y n a m i c M o d e l 44 3 . 2 . 1 . 1 A s s u m p t i o n s 44 3 . 2 . 1 . 2 B a s i c E q u a t i o n s 45 3 . 2 . 2 BOD/DO M o d e l 50 3 . 2 . 2 . 1 A s s u m p t i o n s 50 3 . 2 . 2 . 2 B a s i c E q u a t i o n s 52 3 . 2 . 2 . 3 BOD M o d e l 55 3 . 2 . 2 . 4 DO M o d e l 56 3 . 2 . 3 C o u p l i n g t h e Water Q u a l i t y M o d e l 58 3 . 2 . 4 I n i t i a l and B o u n d a r y C o n d i t i o n s 58 3 .3 N o n l i n e a r Programming M o d e l 61 3 . 3 . 1 S o l u t i o n M e t h o d o l o g y 61 3 . 2 . 2 N o n l i n e a r Programming A l g o r i t h m s 62 3 . 3 . 3 H o o k e - J e e v e s Method 66 3 . 4 Computer R e s o u r c e s and Code 69 4. F I E L D INFORMATION 70 4 . 1 S i t e D e s c r i p t i o n 4 . 1 . 1 C l i m a t e 4 . 1 . 2 G e o l o g y 70 70 72 4 . 1 . 3 R i v e r D i s c h a r g e s , T e m p e r a t u r e and I r r i g a t i 4 . 1 . 4 T i d e s 4 . 1 . 5 L a n d Use and P r i n c i p l e J u r i s d i c t i o n s 4 .2 D a t a S o u r c e s MODEL CALIBRATION AND VALIDATION 5 .1 H y d r o d y n a m i c M o d e l 5.2 BOD/DO M o d e l 5 . 2 . 1 V a l i d a t i o n i n t h e L i t e r a t u r e 5 . 2 . 2 V a l i d a t i o n A g a i n s t A n a l y t i c S o l u t i o n s 5 . 2 . 3 V a l i d a t i o n U s i n g A l t e r n a t i v e R i v e r D a t a 5 . 2 . 4 Summary MODEL A P P L I C A T I O N 6.1 BOD/DO S i m u l a t i o n F a c t o r s 6 . 1 . 1 D i s c h a r g e 6 . 1 . 2 T i d e s 6 . 1 . 3 U n i t BOD L o a d i n g s 6 . 1 . 4 P h o t o s y n t h e s i s 6 . 1 . 5 DO R e g u l a t o r y L i m i t 6 . 1 . 6 M o d e l P a r a m e t e r s U s e d 6.2 T e s t o f U n i m o d a l A s s u m p t i o n 6.3 R e s u l t s and S e n s i t i v i t y A n a l y s i s 6 . 3 . 1 D i s c h a r g e L e v e l s - v i -6 . 3 . 2 DO R e g u l a t o r y L i m i t 6 . 3 . 3 C o m p l i a n c e L e v e l s 6 . 3 . 4 M o d e l P a r a m e t e r s 6 . 3 . 5 O x y g e n S o u r c e s 6 . 3 . 6 O t h e r P o l l u t i o n Abatement P o l i c i e s 7. GENERAL DISCUSSION 7.1 H y d r o d y n a m i c M o d e l 7.2 BOD/DO M o d e l 7 .3 U n c e r t a i n t y o f M e a s u r e d BOD 7.4 D a t a G a t h e r i n g 7 .5 NLP F o r m u l a t i o n 7 . 6 NLP A l g o r i t h m 7 .7 Computer R e s o u r c e s 7 .8 U t i l i t y o f M e t h o d o l o g y 8. CONCLUSIONS REFERENCES - v i i -LIST OF TABLES I Parameters Used for Analytical Validation 94 II Slug Load Comparison 95 III Steady State Comparison 95 IV Sensitivity Ranges for Parameters Tested 111 - v i i i -L I S T OF FIGURES 2 . 1 Components o f t h e Oxygen Sag Curve 8 2 .2 D e t e r m i n i s t i c and Monte C a r l o S i m u l a t i o n 17 2 .3 D e f i n i t i o n o f t h e Sys tem a n d A s s o c i a t e d V a r i a b l e s 18 2 .4 T e m p o r a l and S p a c i a l V a r i a t i o n i n Do L e v e l s 23 2 . 5 T y p i c a l O p t i m a l A e r a t i o n C a p a c i t y A l l o c a t i o n P o l i c y 30 3 .1 S k e t c h o f N o t a t i o n U s e d —47 3 .2 M o d i f i e d Computer Scheme 49 3 .3 F l o w C h a r t o f t h e H y d r o d y n a m i c M o d e l 51 3 . 4 C o n v e c t i o n o f S l u g L o a d 54 3 . 5 F l o w C h a r t o f t h e BOD M u l t i - S t e p P r e c e d u r e 57 3 . 6 F l o w C h a r t o f t h e C o u p l e d M o d e l s 59 3 .7 BOD/DO M o d e l S c h e m a t i c o f R i v e r Reaches 60 3 .8 F l o w C h a r t o f S o l u t i o n M e t h o d o l o g y 63 3 . 8 a ) I n i t i a l i z a t i o n 64 3 . 8 b ) S i m u l a t e BOD C o n c e n t r a t i o n 64 3 . 8 c ) Compute T o t a l BOD C o n c e n t r a t i o n 65 3 . 8 d ) S i m u l a t e DO C o n c e n t r a t i o n 65 3 . 8 e ) P r o p o r t i o n o f C o m p l i a n c e 65 3 .9 Hooke a n d J e e v e s L o g i c D i a g r a m 68 4 .1 L o c a t i o n P l a n s h o w i n g S u r r e y M u n i c i p a l i t y 71 4 .3 P h o t o g r a p h s o f t h e N i c o m e k l R i v e r 75 - i x -4 .4 N i c o m e k l R i v e r P r o f i l e 76 4 . 5 a ) P l a n o f t h e N i c o m e k l - S e r p e n t i n e F l o o d p l a i n 77 4 . 5 b i ) G r a p h o f C r o s s - S e c t i o n a l A r e a v e r s u s Water E l e v a t i o n 78 4 . 5 b i i ) G r a p h o f P e r i m e t e r v e r s u s Water E l e v a t i o n 78 4 . 5 b i i i ) G r a p h o f R i v e r T o p W i d t h v e r s u s D e p t h o f F l o w 79 4 . 6 N i c o m e k l R i v e r D i s c h a r g e s b e l o w M u r r a y C r e e k 80 4 .7 F r e q u e n c y of L o g - N o r m a l i z e d N i c o m e k l D i s c h a r g e s 81 4 .8 S e r p e n t i n e R i v e r S t ream T e m p e r a t u r e s 82 4 .9 T i d a l Gauge H e i g h t a t T s w a s s e n , A u g u s t 1977 84 4 .10 L a n d Use P l a n o f S u r r e y 85 5 .2 V a l i d a t i o n o f t h e H y d r o d y n a m i c M o d e l 91 5 .3 DO P r o f i l e s i n Response t o D i f f e r e n t P a r a m e t e r s 93 5 .4 P04-P C o n c e n t r a t i o n s , O b s e r v e d and S i m u l a t e d 93 5 . 5 a ) D i s s o l v e d O x y g e n a t t h e D e l a w a r e M e m o r i a l B r i d g e 98 5 . 5 b ) D i s s o l v e d O x y g e n a t T o r r e s d a l e , P a . 99 5 . 5 c ) D i s s o l v e d O x y g e n a t B r i s t o l , P a . 100 6.1 F r e q u e n c y o f S t o r m I n t e r a r r i v a l T i m e , N i c o m e k l R i v e r 103 6 .1a B0D/D0 M o d e l S c h e m a t i c o f R i v e r Reaches 106 6.2 G r a p h i c a l R e p r e s e n t a t i o n o f a L i m i t e d G r i d S e a r c h 110 6.3 O p t i m a l BOD L o a d i n g s a t 5 c f s D i s c h a r g e ^ 112 6.4 O p t i m a l BOD L o a d i n g s a t 25 c f s D i s c h a r g e 112 6 .5 O p t i m a l BOD L o a d i n g s a t 3 7 . 5 c f s D i s c h a r g e 112 6.6 O p t i m a l BOD L o a d i n g s a t D = 0 . 5 m i 2 p e r day 117 6.8 O p t i m a l BOD L o a d i n g s a t k i = 0 . 0 5 p e r day 117 6 .9 O p t i m a l BOD L o a d i n g s a t 25% I n c r e a s e i n Oxygen S o u r c e 117 - x -BISMILLAHIR RAHMANIR RAHIM ACKNOWLEDGEMENTS The author extends his appreciation to Dr. M. Puterman, Dr. D. Uyeno, and Dr. W. Caselton for t h e i r guidance and patience during the course of t h i s 'long-distance' t h e s i s . Thanks are also due to Dr. S.O. Russell for his review. Special thanks go the Edmonton Pol i c e Department for the a v a i l a b i l t y of the microcomputers and word processing f a c i l i t i e s for the purposes of t h i s thesis study. This thesis i s dedicated to Zenobia for the long hours spent caring for A l i y a and Khalid during my absences. - 1 -INTRODUCTION As multiple uses of water resources within a r i v e r basin become more competitive, the management of the resource increases i n complexity often aggravated by i n s u f f i c i e n t and uncertain f i e l d information and inadequate control mechanisms. Another order of complexity l i e s i n the dynamic nature of the system being managed. The p h y s i c a l , chemical and b i o l o g i c a l processes need to be well understood i f decisions made at various control l e v e l s and i n various j u r i s d i c t i o n s are to be made with increased confidence. The p o l l u t i o n of downstream r i v e r and groundwater by upstream a g r i c u l t u r a l enterprises has been a well documented, increasing phenomena i n i n t e n s i v e l y developed watersheds. (Tubbs and Haith, 1981) The r i v e r basin manager requires timely an accurate information on the state of the system (Min. of Env., 1984). This information i s often obtained from data gathering networks but increasingly from the judicious use of well validated mathematical models which predict the response of the system to management control acting as amplifiers of f i e l d data. The Nicomekl-Serpentine River Basin i n Surrey, B.C. i s prim a r i l y used for a g r i c u l t u r e which i s a s i g n i f i c a n t contributor of point ans non-point organic p o l l u t i o n . The impact of these animal wastes i n greatest on the f i s h e r y resource which also uses the r i v e r f o r spawning (Moore, 1985a). This t i d a l watershed i s a complex ecosystem where various water q u a l i t y parameters could t y p i f y i t s various s t a t e s . The choice of these state variables i s a n o n - t r i v i a l task e s p e c i a l l y with the multitude of pollutant constituents. - 2 -various states and to develop a mathematical model. The choice of state variables i s a n o n - t r i v i a l task e s p e c i a l l y with multiple pollutants ( f e r t i l i z e r , herbicides and animal wastes). A d d i t i o n a l l y , optimal water q u a l i t y controls need to be based on the e x i s t i n g j u r i s d i c t i o n a l mandates i n the region for maximum effectiveness and enforcement. The problem addressed i n t h i s study i s viewed from a management perpective of determining the to optimal organic loading of a t i d a l r i v e r given i t s capacity for d i l u t i o n and transport of p o l l u t a n t s . S p e c i f i c a l l y , the objective would be to maximize the pollutant loadings along the r i v e r subject to water q u a l i t y regulations. For t h i s study, the basic parameters chosen to describe the state of the system were the organic biochemical oxygen demand (BOD*) and the r e s u l t i n g dissolved oxygen (DO*) concentrations i n the r i v e r due to waste loadings. This reduces the dimensions of the problem to two manageable variables f o r the puposes of t h i s study. A more extensive e f f o r t based on adequate f i e l d data would be required to address the multidimensional nature of the major p o l l u t i o n processes involved. The use of mathematical models i n conjunction with s p e c i f i c management function models (as opposed to research models) i s gaining greater j u s t i f i c a t i o n i n the l i t e r a t u r e (Beck, 1985), and i s the approach adopted here. * Throughout t h i s study, these two parameters w i l l be referred to as BOD and DO res p e c t i v e l y . - 3 -The specific objectives of this study were: 1) A survey of the literature in the areas of water quality modeling for management purposes; 2) Develop a water quality model to account for the major dynamic processes of hydrodynamics and BOD/DO pollution in a t i d a l river; 3) Develop a programming model for the puposes of identifying the optimal BOD loadings subject to regulatory constraints for the purpose of management policy formulation. The literature was surveyed in chapter 2 describing the development and application of models to the BOD/DO water quality problem. Chapter 3 developed the BOD/DO and programming model used to generate optimal policies in this study. The BOD/DO model is validated in chapter 5 and applied in chapter 6, while chapter provides background site information. Discussion and conclusions follow in chapters 7 and 8 including suggestions for further work. 2. LITERATURE SURVEY 2.1 INTRODUCTION The a p p l i c a t i o n of systems engineering techniques to the solution of q u a l i t y problems of water bodies i s an ongoing area of research. Reviews of the f i e l d (for example the Journal of the Water P o l l u t i o n Control Federation annual review) document the large e f f o r t i n academia, government and industry i n the development of suitable techniques and methodologies i n r e s o l v i n g user c o n f l i c t s and improving the management of water bodies. This growing p o r t f o l i o of techniques can be used to a s s i s t a manager r a t i o n a l i z e suitable control s t r a t e g i e s , at a minimum cost of information and c o n t r o l , based on knowledge of p h y s i c a l processes within the water resource system. Models of complex systems are both d i f f i c u l t to b u i l d and to v a l i d a t e but are necessary i f a system i s to be m u l t i - v a r i a b l y c o n t r o l l e d . In other words the model should have the 'requisite v a r i e t y ' (Beer, 1966). However, i n most cases f i e l d data are lacking and only the use of s i m p l i f i e d e c o l o g i c a l models, such as a BOD/DO model, can be j u s t i f i e d using parameter estimation and cautious f i e l d a p p l i c a t i o n . The need also e x i s t s to e x p l i c i t l y model the stochastic nature of the water q u a l i t y problem in sampling, management and enforcement of regulations. Ward and L o f t i s (1983) state: "Water q u a l i t y involves a mix of deterministic and stochastic concepts, and i t can be e f f e c t i v e l y managed only when both components are properly balanced". This chapter surveys methods ava i l a b l e i n the l i t e r a t u r e i n the optimal management and control of r i v e r water q u a l i t y . - 5 -2.2 WATER QUALITY MODELS Beck (1983) proposed the following c l a s s i f i c a t i o n of model types: i ) research versus management models; where the former provides ind i c a t o r s for f r u i t f u l d i r e c t i o n s for i n v e s t i g a t i o n while the l a t t e r are for s t r a t e g i c and t a c t i c a l water q u a l i t y management. i i ) d i s t r i b u t e d versus lumped models; where the former attempts to model processes of i n d i v i d u a l constituents i n a water q u a l i t y system while the l a t t e r aggregates these v a r i a t i o n s . i i i ) non-linear versus l i n e a r models; where i n the former an increase of an input i s r e f l e c t e d i n a non-linear model response whereas the l a t t e r predicts a l i n e a r response. iv) stochastic versus deterministic models; where the former accounts for the random inputs into the system which are ignored by the l a t t e r . v) dynamic versus steady-state models; where the former models the time dimension of the system while the l a t t e r does not. vi ) i n t e r n a l l y d e s c r i p t i v e versus 'black box' models; where the former deals with a p r i o r i information and deductive reasoning based on knowledge of the processes involved while the l a t t e r i s p o s t e r i o r i and inductive based on observation. - 6 -Within t h i s c l a s s i f i c a t i o n we s h a l l survey lumped parameter (in one dimension), nonlinear and l i n e a r , stochastic and det e r m i n i s t i c , dynamic and steady-state, and i n t e r n a l l y d escriptive models. Most f i e l d a pplications of BOD/DO models are deterministic, steady-state and one-dimensional while unsteady three-dimensional stochastic e c o l o g i c a l models tend to be research oriented. 2.2.1 Steady-State Models The most common BOD/DO model i s that of Streeter-Phelps which assumes that: a) the BOD decay rate i s proportional to BOD concentration; b) the deoxygenation and BOD decay rate are equal; c) the reoxygenation rate i s proportional to the oxygen d e f i c i t ; d) BOD and DO are s u f f i c i e n t to describe the biochemical process| r e s u l t i n g i n the following non-linear model: db + vtb = -k..b (2.1a) b t 61 be + v6c = - k ^ + k 2 ( C g - c) bt 61 (2.1b) with boundary conditions: b ( 0 , t ) = b Q ( t ) and c(o,t) = c Q ( t ) , and i n i t i a l conditions: b ( l , 0 ) = b ^ ( l ) and c ( l , 0 ) = c ^ ( l ) , THe boundary conditions can vary i n any manner from impulse functions through sinusoidal functions to steady-state values. This i s also true for the i n i t i a l conditions along the stretch of the r i v e r considered. In equations 2.1: b ( l , t ) - the (ultimate) BOD - 3 concentration i n the r i v e r [M.L ] c ( l , t ) - the DO concentration [M.L~3] t - i s the time dimension 1 - i s the space dimension V - i s the average stream v e l o c i t y [L.T"1] k1 - the deoxygenation c o e f f i c i e n t [T~1] k2 - the reaeration c o e f f i c i e n t [T"1] C s - saturation DO concentration [M.L~3] three parameters are temperature dependent as k T= k20 e < The BOD i s normally evaluated as B0D n» where n i s usually taken as f i v e days and BOD = (1- e~ K L n)B0D , with KT the BOD 1 n ' u l t L decay c o e f f i c i e n t derived from standard laboratory t e s t s . Solution of equations 2.1 with steady-state assumptions (b (t)=b ,c (t)=c / o o o o bb/dt = 0, 6c/6t = 0) y i e l d s the c l a s s i c a l Streeter-Phelps equation: b(T) = b Q e " k - | T (2.2a) c(T) = C -(C - c 0 ) e " k 2 T + b ( e " k l T - e ~ k 2 T ) K1 K2 (2.2b) where dl/dT = v and b Q r e s u l t s from both carbonoceous and nitrogeneous BOD loadings and k^ i s a composite deoxygenation c o e f f i c i e n t . (T) and A l t e r n a t i v e l y , b(T) could be computed as b carb , (T) re s p e c t i v e l y . Equation 2.2 describes the oxygen-sag curve n i t r shown i n figure 2.1. Dobbins (1964) proposed modifications based on the following, more e x p l i c i t , assumptions: a) the stream flow i s steady and uniform; - 8 -o ~ 0 ieficit 0 2 u ^ger (mg/lii 4 (xo ) (mg/lii 6 Ne 8 D a y s after ef f luent a d d e d to s t r e a m D i s t a n c e b e l o w o u t f a l l *-Figure 2.1 C o m p o n e n t s o f the o x y g e n sag curve : (a) D e o x y g e n a t i o n , the rate o f w h i c h is p r o p o r t i o n a l to -the c o n c e n t r a t i o n o f r e m a i n i n g B O D i n the water , t h e r e f o r e d e c l i n e s as the o x i d i z a b l e m a t e r i a l is c o n -s u m e d , (b) R e a e r a t i o n , the rate o f w h i c h is p r o p o r t i o n a l to the o x y g e n de f ic i t at the t ime , (c) N e t o x y g e n def ic i t o b t a i n e d b y s u b t r a c t i n g reae ra t i on f r o m d e o x y g e n a t i o n . T h e f o r m o f the cu rve leads to its b e i n g c a l l e d the " o x y g e n s a g . " (Dunne and Leopold. 1978) - 9 -b) the process for the stret c h or r i v e r as a whole i s steady-state, the condition at each cross-section being unchanged with time; c) the removal of BOD by both b a c t e r i a l oxidation and sedimentation and/or adsorption are f i r s t - o r d e r reactions, the rates of removal at any cross-section being proportional to the ammount present; d) the removal of oxygen by benthal demand and plant r e s p i r a t i o n , the addition of oxygen by photosysnthesis, and the addition of BOD from benthal layers or l o c a l runoff are c a l l uniform along the str e t c h ; e) the BOD and oxygen are uniformly d i s t r i b u t e d over each cross-section thus permitting the equations to be written i n the usual one-dimensional form. These assumptions suggest the following non-linear steady-state model: D Td 2b - vdb + (k.+k-Jb + b = 0 (2.3a) dx dx D Ld 2^ + k 2 ( C g - c ) - vdc - k ^ - D B = 0 (2.3b) dx dx where: 2 DT - the lo n g i t u d i n a l dispersion c o e f f i c i e n t [L .T] k^ - c o e f f i c i e n t of BOD removal by - 1 sedimentation/adsorption [T ] D„ - net rate of removal of oxygen by benthal demand and photosynthesis [M.T b - rate of addition of oxygen along - 10 -the stre t c h by a r t i f i c i a l aeration [M.T~^] The solution of equation 2.3 , with n e g l i g i b l e l o n g i t u d i n a l dispersion, i s given as: b~ b©6 ^'+^+b_ ( l ~ C^*^^) ( 2 , 4 a) and: Q-C - k, ( b 0 - k ) r e - < k - + l ^ f c . e ] (2.4b) I 2_ + k,b_ V i - e ) where b and c are i n i t i a l BOD and DO concentrations, o o A v a r i e t y of other one-dimensional models have been formulated depending on the assumptions made. For example Shastry et a l (1973) proposed a non-linear decay equation as _rl? * - k p b e 4. k»,b •* i ^ j , and e_ = __££• + k_ccs-c) + D & with as a saturation constant. These equations combined never predict negative DO values as does the Streeter-Phelps model. This class of steady-state models are r e l a t i v e l y easy to combine, d i r e c t l y or i n d i r e c t l y , with optimization procedures such as mathematical programming or optimal control and to apply to problems of resource a l l o c a t i o n , treatment or regulation where BOD and/or DO loadings could be co n t r o l l e d . The s e n s i t i v i t y of these models and the r e s u l t i n g control p o l i c i e s i s dependent on the model parameters chosen. Wen et a l (1982) derived s e n s i t i v i t y functions for BOD and DO for a Dobbin's type model to decay rate and dispersion parameters i n non-tidal and t i d a l r i v e r s where the c o e f f i c i e n t of l o n g i t u d i n a l dispersion was zero and non-zero re s p e c t i v e l y . The a p p l i c a t i o n of steady-state type models i s extensively documented i n the l i t e r a t u r e due to t h e i r r e l a t i v e s i m p l i c i t y , and hence greater p o t e n t i a l for mis-use, since the correct determination of parameter values i s a n o n - t r i v i a l task. Rickert et a l (1976) applied a steady-state model (incorporating nitrogenous BOD demand) to the Willamette River, Oregon to t e s t the impact of management al t e r n a t i v e s such as BOD loading a l l o c a t i o n s , ammonia loadings, low flow augmentation and the reduction of benthos on r i v e r DO. Ramm (1976) determined the impact of i n d u s t r i a l and domestic wastes of the DO p r o f i l e s i n the navigation channel of the Cuyahoga River. They used a compartment model (where mass balances were conserved over r i v e r reaches as opposed to elemental sections) f o r extensive s e n s i t i v i t y analysis and presented a t r a n s f e r matrix for use as a loading decision table. Putz et a l (1983) used a steady-state transverse mixing model to predict the two-dimensional dispersion and decay of f e c a l coliforms from a proposed lagoon for Fort Smith on the Slave River, NWT. Their model, va l i d a t e d with steady dye tracer studies, predicted higher than allowed coliform concentrations during winter. Kingscott (1976) used a steady-state model to assess the cost effectiveness of waste load a l l o c a t i o n s to various treatment unit processes by p r e d i c t i n g the r e s u l t i n g r i v e r DO due to stochastic (Monte Carlo generated) BOD loadings. Hung et a l (1976) presented a modified Streeter-Phelps model used to e s t a b l i s h waste load a l l o c a t i o n s for Indiana r i v e r s by combination with optimization models and p r o b a b a l i s t i c - 12 -methods. Yearsley (1976) also discussed Monte Carlo techniques used to generate DO p r o f i l e s for the Milner reach of the Snake River, Idaho. The inherent v a r i a b i l i t y i n model parameters (recognized and accounted for i n the l a s t three references above) has prompted the development of a n a l y t i c stochastic•steady-state models where these parameters are considered to be stochastic v a r i a b l e s . Thayer and Krutchkoff (1967) formulated a d i f f e r e n t i a l difference equation describing the p r o b a b a l i s t i c dynamics of BOD and DO in t e r a c t i o n s which they solved using moment generation functions i n deriving j o i n t p r o b a b i l i t y d e n s i t i e s . They validated t h e i r model by suc c e s f u l l y matching DO variances for the Sacramento River with t h e i r predicted values. Padgett and Durham (1976) solved a p r o b a b i l i s t i c version of Dobbins' model as a set of "random d i f f e r e n t i a l equations with random c o e f f i c i e n t s , random inhomogeneous terms, and random i n i t i a l conditions' which they solved using the mean square theory approach proposed by Soong (1973). Further work by Padgett et a l (1976), B e l l and Papadopolous (1978), Padgett and Papadopolous (1978) and Damaskos and Papadopolous (1983) has led to a formulation of a general stochastic BOD/DO model. These models y i e l d stochastic steady-state BOD/DO p r o f i l e s f o r various parameter j o i n t d i s t r i b u t i o n a l assumptions such as normality and independence with stochastic i n i t i a l and boundary conditions. These models however, do not appear to have been f i e l d tested. 2.2.2 Dynamic Models Deterministic water q u a l i t y models which describe the ph y s i c a l , chemical and b i o l o g i c a l processes found i n water bodies are widely used i n the assessment of impacts of pollutant loadings. However, the app l i c a t i o n of such models can, as mentioned previously, be severely - 13 -l i m i t e d by the lack of f i e l d data, the lack of knowledge of the underlying processes, and the lack of advanced modeling techniques. (Mausberger, 1983). The general mass balance of a constituent i n a segment of water i s Rate change of + Rate of + Rate of = Net Rate of Mass Density Convection D i f f u s i o n Addition In one-dimension t h i s reduces to the advection-diffusion equation: MAO _ W A E k ^ . - M Q < 0 + A G (2.5) 2 where A - cr o s s - s e c t i o n a l area [L ] _3 c - concentration of the constituent [M.L ] 3 — 1 Q - r i v e r discharge [L .T ] 2 -1 E - e f f e c t i v e dispersion c o e f f i c i e n t [L .T ] _ i G - source-sink term [M.T ] t - time [TJ x - distance [L] For an organic pollutant (ultimate BOD) to equation takes the form: M A b ) _ b_(AB_D\ _ j _ _ ( Q b ) - ( K ' ^ ) A b + A b a 'bt b y ^ &X (2.6) _3 where b - the ultimate f i r s t - s t a g e BOD [M.L ] k.j - decay c o e f f i c i e n t for BOD removal . [T ^ ] k^ - decay c o e f f i c i e n t for BOD removal by sedimentation and adsorption [T ^ ] L - rate of BOD addition along the stream [M.T 1] cl - 14 -For dissolved oxygen (DO) the equation i s : 6 (Ac) K b £>x^ 5 * / b * (2.7) where: c - DO concentration [M.L - 3] C g - saturated DO concentration [M.L - 3] reaeration c o e f f i c i e n t [ T _ 1 ] D„ - net rate of DO addition due to photosynthesis, r e s p i r a t i o n and benthic demand [M.L~ 3.T~ 1] (Dresnack and Dobbins, 1968) Using v a r i a t i o n s of the above mathematical models, many studies reported i n the l i t e r a t u r e have assumed steady-state hydrodynamics (Q=QQ) while pollutant loadings vary. P h e i f f e r et a l (1976) apply a dynamic BOD model coupled with a quasi-dynamic hydraulic model, which generated stepped streamflow hydrographs, to the Patuxent River Basin to determine the e f f e c t of secondary treatment plants on stream DO. Pence et a l (1968) apply a dynamic BOD/DO model with steady r i v e r flows to the Delaware Estuary to evaluate control schemes to a l l e v i a t e low DO conditions due to municipal, i n d u s t r i a l , t r i b u t a r y and stormwater overflow BOD loadings. Weatherbe (1976) used a s i m i l a r model i n es t a b l i s h i n g water q u a l i t y guidelines for the Thames River, England using h i s t o r i c a l data to generate d a i l y streamflow by 'stochastic' techniques. the r i v e r . Assuming that the i n t e r a c t i o n between water q u a l i t y e f f e c t s and water quantity processes are independent, namely, that constituent Other studies, however, also accounted for the hydrodynamics of - 15 -concentrations do not a f f e c t the density of r i v e r water, the model can be decoupled with the hydrodynamics described by the equations of continuity and momentum. Continuity: (2.8) 0^ + UdU = - c x d h - a 10| U where U - stream v e l o c i t y i n d i r e c t i o n x [L.T~^] g - acc e l e r a t i o n due to gravity (9.81 m.sec ) c - Chezy c o e f f i c i e n t h - flow l e v e l with respect to datum [L] d - depth of flow [L] Momentum: (2.9) a (AO) 4 b 3h = o dx "St where b - top r i v e r width [L] (Dronkers, 1969). Tucci and Chen (1981) used the Saint Venant formulation of the equations of continuity and momentum above, coupled with a one-dimensional transport equation i n modeling the e s t u a r i a l r i v e r network of the Jacui Delta, B r a z i l for the purpose of evaluating management options for BOD loadings i n the four t r i b u t a r y r i v e r s . Johanson et a l (1976) applied a noded coupled unsteady model to the harbour and lower Chahalis River, Washington. They conducted s e n s i t i v i t y analyses of on parameter values while evaluating the DO impact of i n d u s t r i a l discharges. Sargent (1976) used plug-flow and completely-mixed coupled models (with zero and non-zero dispersion respectively) i n evaluating the impact of water treatment on the Buffalo River. Bansal (1976) applied dynamic and steady-state models to the Big Blue River i n Nebraska concluding that the one-dimensional dispersion model was the most s a t i s f a c t o r y i n - 16 -v e r i f i c a t i o n with f i e l d data. Beckers et a l (1976) applied coupled models to predict the concentration of phosphorus, the 'nitrogens', coliforms, BOD, DO, chlorophyll and s a l i n i t y to channels of the Pawtuxet River of Rhode Island where f i e l d data was s u f f i c i e n t , concluding that the a p p l i c a t i o n of unsteady models was preferred i n a p p l i c a t i o n to dynamic processes found i n estuaries and stormwater runoff. The parameters i n a dynamic model are, as mentioned previously, stochastic i n nature either due to the lack of understanding of the phenomena involved or to inherent v a r i a b i l i t y of the natural processes. Since equations 2.6 and 2.9 are solved numerically, Monte Carlo simulation methods are commonly used often with Markovian generation techniques. Ford et a l (1980), for example, used such methods i n developing 95% confidence envelopes for the p r e d i c t i o n of a l g a l concentrations i n the proposed Twin Valley lake using the U.S. Corps of Engineers CE-QUAL-RI model as shown i n figure 2.2. C o e f f i c i e n t values were independently generated from a v a r i e t y of p r o b a b i l i t y d i s t r i b u t i o n s using recorded mean and variance data. Similar methods were used i n a study of the water q u a l i t y impacts of the proposed B.C. Hydro dams on the Stikine and L i a r d Rivers of northern B r i t i s h Columbia (Schultz, 1983). Scholl and Wycoff (1981) also applied Markov models i n generating c o e f f i c i e n t values for t h e i r Continuous Stormwater P o l l u t i o n Simulation System i n t h e i r analysis of the impact of the S p r i n g f i e l d advanced water treatment plant on the James River i n Missouri and the i d e n t i f i c a t i o n of the p o l l u t a n t sources responsible for f i s h k i l l s . - 17 -is 2 <J> 9 6 5 M 3 / S 4 8 M 3 / S I29 159 I89 2 I 9 2 J U L I A N DAY 19 75! a. DETERMINISTIC SIMULATION 2 7 9 3 0 9 1 0 8 ///;/;.e& M 3 / S S S ^ S 4 8 M 3 / S 1 9 8 2 2 8 2 5 8 J U L I A N D A Y 1 9 7 5 b. MONTE CARLO SIMULATION 2 8 8 3 1 8 F I G . 2.2 C o m p a r i s o n o f A l g a e C o n c e n t r a t i o n s f o r D e t e r m i n i s t i c S i m u l a -t i o n a n d M o n t e C a r l o G e n e r a t e d 95 P e r c e n t C o n f i d e n c e I n t e r v a l f o r T w o R e l e a s e R a t e s . (Ford et al, 1980) - 18 -2.2.4 Estimation Techniques Parameter estimation techniques account for the v a r i a b i l i t y i n parameters of c o r r e c t l y i d e n t i f i e d and v e r i f i e d water q u a l i t y models (Beck, 1985). The system can be viewed as i n figure 2.3 where uncertainty i s found i n both the processes observed and i n the observation process, Unmeasured d is turbances Measurement errors Measured input d is turbances Process States x Parameters a Vy Measured o u t p u t s FIGURE 23Definition of the system and associated variables. (Beck, 1985) where the v a r i a t i o n s with times of the system state vector are given by dx(t) = f [ x ( t ) , u ( t ) , a ( t ) , e(t) ] (2.10) dt with sampled discrete-time output observation y/tj^) = h [ x ( t k ) , a ( t k ) + n ( t k ) ] (2.11) where x = n-dimensional vector of state variables u = m-dimensional vector of measured input disturbances _y_ = 1-dimensional vector of ( d i s c r e t e l y sampled) measured output variables a = p-dimensional vector of model parameters j 2 = s-dimensional vector of random unmeasured (unknown) input disturbances n = 1-dimensional vector of random output measurements - 19 -and f and h are nonlinear, vector-valued functions; t i s the independent time dimension and t ^ i s the k'th discre t e sampling instant i n time. (Beck, 1985). Two major estimation e f f o r t s a r i s e : state estimation and parameter estimation whereby optimal estimates of vectors _x and a_ are made using, for example, maximum liklehood, least squares, Markov and Bayesian estimators. However, a s i g n i f i c a n t problem that arises i s that of ' i d e n t i f i a b i l i t y ' , whereby a model i s chosen which 'best' simulates the output of the system (Rinaldi et a l , 1979). Parameter estimation techniques are used both ' o f f - l i n e ' and i n 'real-time' to estimate steady as well as varying parameters. Tamura and Kawaguchi (1980) presented a real-time parameter estimation procedure for BOD parameters i n four reaches of the Yomo River, Japan. The procedure was based on a method f i t t i n g n_-dimensional time series data to an autoregressive model of f i n i t e order, using a steady-state model. Olsen (1980) applied the techniques of dynamic modeling, real-time p r e d i c t i o n , estimation and control to the Kappala waste-waster treatment plant, Sweden, concluding that the p r e d i c t i o n and estimation methods were probably the most p r o f i t a b l e areas for the ap p l i c a t i o n of control theory to waste-water treatment systems. Rinaldi et a l (1979) reviewed three major estimation a p p l i c a t i o n s : a) parameter estimation of a modified Streeter-Phelps model using a d i r e c t search method for a nonlinear programming formulation of minimized squared deviations; b) state and parameter estimation by the q u a s i - l i n e a r i z a t i o n technique for solving nonlinear multi-point boundary value problems. - 20 -Squared deviations, incorporating measurement noise, were minimized using maximum l i k e l i h o o d estimators for the Rhine River; c) Kalman f i l t e r i n g techniques were applied to state and parameter estimation for the Yomo River, based on a steady-state, d i s t r i b u t e d - l a g model and compared with other f i l t e r s i n terms of e f f i c i e n c y . These researchers concluded that a majority of water q u a l i t y studies apply a p r i o r i developed models to water q u a l i t y systems with the purpose of showing t h e i r i n v a l i d i t y , rather than using p o s t e r i o r i techniques of estimation to i d e n t i f y the appropriate model. On t h i s same theme Beck (1983, 1985) noted a trend towards po s t e r i o r analyses as the f i e l d of water q u a l i t y modeling matures. An example of such an analysis i s the work of Shastry et a l (1973) who used nonlinear estimation techniques to s e l e c t the most appropriate steady-state model (of three considered) i n the simulation of DO i n the Sacramento River. Applying least squares and maximum liklehood methods, they estimated both model parameters and i n i t i a l conditions for the r i v e r and noted that the unimodality of the nonlinear response surface was d i f f i c u l t to determine a p r i o r i , e s p e c i a l l y for high dimensioned parameter vectors. 2.3 SYSTEMS OPTIMIZATION TECHNIQUES The a p p l i c a t i o n of optimization methods to the control of water p o l l u t i o n has s t e a d i l y developed since the early applications of l i n e a r programming i n the 1960's. These methods can be generally c l a s s i f i e d as: - 21 -1) Linear Programming (LP) 2) Nonlinear Programming (NLP) 3) Dynamic Programming (DP) 4) Stochastic Programming (SP) 5) Optimal Control (OC) 2.3.1 Linear Programming Formulations A t y p i c a l l i n e a r program i s (max) C'X (min) such that AX <_ b and X > 0 where the objective function optimizes, for example, treatment costs or pollutant loadings. The constraints describe budget r e s t r i c t i o n s , capacity l i m i t a t i o n s or pollutant discharge regulations. This formulation i s d i s c r e t i z e d over space where the constraints could apply to a f i n i t e set of monitoring s t a t i o n s . The cost function may also require to be l i n e a r i z e d r e s u l t i n g i n large LP models. Rinaldi et a l (1980) suggest a decomposition technique to reduce the LP to a multistage problem with a reduction i n the number of constraints by i d e n t i f y i n g c r i t i c a l monitoring points and the use of l i n e a r approximations to the steady-state DO sag l e v e l s as, for example, proposed by Arbabi et a l (1974). This approximation r e s u l t s i n a DO based on the Streeter-Phelps model. Biswas and Reynold (1973) used l i n e a r i z e d Streeter-Phelps equations i n a screening model i n i d e n t i f y i n g p ossible minimum-cost treatment operating points i n t h e i r study of DO regulation of the Saint John River Basin of New Brunswick, Quebec and Maine. They then used a dynamic pollutant model, driven by synthetic r i v e r hydrographs, to predict seasonal DO p r o f i l e s as shown i n f i g u r e ion - 22 -2.4. Loucks (1983) discussed the following problems for solution by l i n e a r programming or mixed integer formulations: a) regional treatment and transport where t o t a l costs are minimized subject to conservation of flow constraints; b) multiple point source waste reduction to meet water q u a l i t y standards where costs of treatment were minimized subject to BOD/DO dynamics modeled ( l i n e a r i z e d s i m i l a r l y to Arbabi et a l , 1974); c) p r e d i c t i o n of a l g a l bloom p o t e n t i a l by maximizing a l g a l production subject to nutrient and l i g h t c onstraints. Due to the steady-state assumptions required for the LP, these approximate formulations are often used as screeening models for further d e t a i l e d simulation models which are then used to explore the optimal sol u t i o n . The dynamics of the r i v e r processes, however, would have been averaged i n a r r i v i n g at t h i s optimum. An a l t e r n a t i v e approach which incorporates the pollutant dynamics into an LP formulation was proposed by Gorelick and Remson (1982) and Gorelick (1982). This methodology i s applicable to any l i n e a r , d i s t r i b u t e d parameter, solute transport water q u a l i t y model. Gorelick (1982) used the deterministic steady-state groundwater flow equation: where ^ i j - the t r a n s m i s s i v i t y tensor H - the hydraulic head W - the volume f l u x per unit area x. - the two horizontal dimensions [ L 2 . T _ 1 ] [L] [M.L _ 2.T - 1] [L] - 24 -and the deterministic advective-dispersive equation: = d_ (Dtj dc") - d.(OvO - cSy (2.13) at dxj. dx* 0 b —3 where C - concentration of dissolved chemical specie [M.L ] ' 2 —1 D^ _. - d i s p e r s i v i t y tensor [L .T ] - average pore water v e l o c i t y i n the d i r e c t i o n of x^ [L.T ^ ] b - saturated aquifer thickness [L] (fi - e f f e c t i v e aquifer porosity C - solute concentration i n a f l u i d source or sink [M.L 3 J t - time [T] These equations describe the movement of a non-conservative pollutant i n a two-dimensional horizontal aquifer. Unit pollutant slug loadings at point sources i n the aquifer were modeled at observation wells as 'break-through' curves which were assumed to occur annually. These curves were d i s c r e t i z e d over a management period of f i v e years to form a unit concentration response matrix to be used i n the following LP formulation: max Z = u_' f_ subject to R f_ <^  c_*v for each observation well, and f > 0 where R i s the concentration response matrix _f i s the vector of unit pollutant loading rates c* i s the vector of water q u a l i t y standards u' i s the unit vector - 25 -The objective was to maximize the t o t a l , equally weighted loading of the aquifer. This formulation was subsequently solved as a dual LP. This methodology was applied by Danskin and Gorelick (1985) to a multi-aquifer groundwater and surface system near Livermore, C a l i f o r n i a where the cost of water a l l o c a t i o n was minimized subject to physical and economic constraints. Without t h i s management model, these researchers state that the problem posed would have been impossible or extremely d i f f i c u l t to solve. A s i m i l a r methodolgy had also been applied by Bishop et a l (1976) to a steady-state r i v e r with multiple p o l l u t a n t s , where the objective was to minimize treatment costs subject to water q u a l i t y and treatment contraints. These researchers coupled a l i n e a r i z e d steady-state simulation model with an integer programming model v i a a c o e f f i c i e n t matrix and also remarked that i t e r a t i v e s o l u t i o n techniques would have been required for non-linear water q u a l i t y models. 2.3.2 Non-Linear Programming Models These programming models can accommodate the non-linearity of water q u a l i t y processes. These models are usually of the form: where f(.) and g(.) can be convex and/or concave functions, defining either convex or non-convex so l u t i o n sets, to be solved by appropriate numerical techniques. R i n a l d i et a l (1980) discussed the following examples: subject to: g(x) <_ b a) the minimization of waste-water treatment and water production costs less r e c r e a t i o n a l and aesthetic benefits subject to dynamics of a d i s c r e t i z e d Dobbins model - 26 -and BOD/DO standards using the sequential unconstrained minimization technique (SUMT); b) the minimization of treatment costs i n an estuary described by a d i s c r e t i z e d Streeter-Phelps model using Rosen's gradient p r o j e c t i o n method for optimization; c) the minimization of quadratic aeration costs with a Streeter-Phelps model, using v a r i a t i o n a l approaches for s o l u t i o n ; d) a non-convex minimization of non-linear costs of treatment and flow augmentation subject to flow conservation and q u a l i t y standards using a d i r c t search method (ensuring, however only a l o c a l optimum). Bayer (1972) obtained minimum treatment plant, dam and r e s e r v o i r costs subject to BOD and DO standards, low flow l i m i t s and plant e f f i c i e n c y ranges. Both f i r s t - o r d e r (carbonaceous) and second-order (nitrogenous) BOD uptake demands were modeled with d i s c r e t i z e d Streeter-Phelps equations applied to the Willamette River, Oregon. Comparisons were made with previously developed LP and DP models for the same r i v e r showing high agreement i n optimal p o l i c i e s generated. Deininger (1972) discussed a minimum-cost regional p o l l u t i o n control system with the optimal l o c a t i o n of treatment plants and a sewer system. The r e s u l t i n g NLP with concave constraints was solved i t e r a t i v e l y by a procedure in v o l v i n g the formulation of a r e l a t e d LP and a ranking extreme point approach to avoid l o c a l minima. This method of s o l u t i o n was found to be appropriate for 'small to medium' sized problems. Taylor and Judy (1972) used a three phase approach i n solving a hypothetical water q u a l i t y management problem by rent allotment c o n t r o l . An NLP, minimizing discharges plus - 27 -agency costs subject to regulatory standards, was used to generate Kuhn-Tucker conditions which described the optimal use of treatment and flow augmentation. These conditions were then used as a basis for an agency-discharger game, with bidding, which was 'played out' when a set of agreements on rents and allotments had been reached. These agreements were then used i n a rent allotment market model to allow decision makers to adopt more e f f i c i e n t waste abatement measures. These examples dealt with steady-state approximations to dynamic systems with reasonably simple models. The d i s c r e t i z e d version of the dynamic model of section 2.2.3 incorporated into a optimization model, though more de s c r i p t i v e , would r e s u l t i n large NLP or LP models i n order to meet s t a b i l i l t y requirements and would require large scale optimization techniques, such as those proposed by Lasdon (1970), f o r s o l u t i o n . Loucks (1983) describes the conversion of the U.S. Environmental Protection Agency's QUAL II model into an economic management model by e x p l o i t i n g f i n i t e - d i f f e r e n c e techniques for the formulation of a multi-constituent optimization. Similar discussions r e l a t e d to lake and reservoir management were also presented. 2.3.3 Dynamic Programming Formulation Due to the sequential multi-stage nature of the r i v e r water q u a l i t y control problem, i t i s amenable to a c l a s s i c a l dynamic programming formulation as: nAX^ = min [C ± (__._) + H (£_ ±{X ± tu. ±))] u± 1 (2.14) where C^Hi^ ^ s s o r a e cost function, such as treatment or pollutant impact on f i s h . - 28 -i s the state value of the system (for example BOD and DO) _u^  i s the control action taken, such as treatment, aeration, pollutant loadings. fj.(x_/U_.) * s ^ e s t a t e a t t n e e n d °f each stage (reach) i f control u^ i s taken with state values x^ which s a t i s f i e s the constraints. For example, the use of the Streeter-Phelps or input-output models. H.(x.) i s the (minimum) cost function at state x. and 1—1 — l stage i . Due to the r e l a t i v e l y complex nature of the BOD/DO models i n the formulation above, a closed form sol u t i o n i s d i f f i c u l t to obtain, suggesting the use of d i s c r e t e , dynamic programming techniques. The major l i m i t a t i o n i n the use of t h i s method i s the well known 'curse of dimensionality' which l i m i t s a p p l i c a b i l i t y on even the most powerful computers. This has led to the development of DP techniques such as d i f f e r e n t i a l DP, discr e t e d i f f e r e n t i a l DP and state incremental DP which, according to Yakowitz (1982), have the p o t e n t i a l to avoid or eliminate the 'curse'. For deterministic problems, forward or backward stage c a l c u l a t i o n methods can be applied, while stochastic problems require the l a t t e r . The scope of research and a p p l i c a t i o n of DP to water resource problems has expanded considerably i n the past decades with a primary focus on rese r v o i r management (Yeh, 1985). Applications to water q u a l i t y problems have also received increased attention. Chang and Yeh (1973) looked at the problem of optimal a l l o c a t i o n of a r t i f i c i a l aeration along - 29 -a p o l l u t e d stream to to meet a dissolved oxygen standard of 5.0 mg/1. They used the discrete DP approach i n a l l o c a t i n g aerator capacity to equidistant aerators at minimum cost, subject to some t o t a l aeration minJ(D,V) = W.,(DL - D t i! capacity. The formulation used was N ) 2 + w 2 v t . 2 (2.15) 1=1 where and W2 are weighting factors and D L i s the regulated DO d e f i c i t (from saturation) and Dt^ i s the actual DO l e v e l , while V t^ are the aerator DO c a p a c i t i e s , subject to 1) the l i m i t a t i o n of t o t a l a v a i l a b l e aeration capacity; N 2) the physical l i m i t a t i o n of DO d e f i c i t and aeration capacity 0 < D t i < c s 0 < V+.. < C i=1 ,. . .. ,N 3) the DO d e f i c i t at the end of every reach (aerated or not) derived from the Streeter-Phelps equations: Using a Lagrange m u l t i p l i e r technique these researchers added the second constraint to the objective functional while constraints 1 and 3 above defined the admissable region for the problem. Using a range 24 of m u l t i p l i e r s for a 24 stage, 11 state formulation, a l l 11 possible states were searched r e s u l t i n g i n control t r a j e c t o r i e s shown i n f i g u r e 2.5. The Lagrange m u l t i p l i e r was interpreted as the cost of a l l o c a t i n g a unit of aeration capacity to the stream with a zero m u l t i p l i e r value at f u l l capacity. - 30 -^ I . O r £ > o 2.0 -3.0 E 5.0 \ \ X-0.1 - \ w,-0.7 w2=0.3 " • * ~ V ~ " " l , , I 2 3 4 t(DAYS) Figure 2.5 Typical Optimal Aeration Capacity Allocation Policy and the Corresponding DO Sag Profile (Chang and Yeh, 1973) - 31 -Rin a l d i et a l (1979) discussed two DP a p p l i c a t i o n s : one being the optimal a l l o c a t i o n of wastewater treatment to the Rhine River, Germany and the other being the optimal a l l o c a t i o n of a r t i f i c i a l instream aeration i n a reach of the same r i v e r . Both applications are formulated as the optimal a l l o c a t i o n of complementary p o l l u t i o n control mechanisms using deterministic, steady-state BOD/DO models to bound the admissable so l u t i o n region. River conditions were established for the most unfavourable water q u a l i t y conditions during which environmental standards of 4.0 mg/1 DO and 20 mg/1 BOD were to be enforced. The 6 4 treatment problem with 6 states was reduced to 6 states by ex p l o i t i n g the sp e c i a l structure of the DO/BOD model and solved using backward recursion. The aerator problem was formulated as follows: a n d S.t. N d t J -oCe) (Ufa with CXO) = C o , (2.16) where C. i s the cost of aeration [$] Q. i s the r i v e r discharge [ L 3 . T - 1 ] z^(l) i s a general organic decomposition function equivalent to BOD i n the Streeter-Phelps model [M.L - 3] w(l) i s the net source of pollutant 6w i s an impulse function N i s the number of aerators u^ i s the aerator capacity L. i s the stream reach length [L] - 32 -The problem as formulated was s i m p l i f i e d by decomposition i n t o independent sub-problems based on segments of the r i v e r . The water q u a l i t y model was solved using a Runge-Kutta in t e g r a t i o n scheme. An i n t e r e s t i n g extension made was where r i v e r flow was f i t t e d with a piece-wise l i n e a r function by distance along the r i v e r . This was an attempt to describe the spa c i a l dynamics of the r i v e r and would, however, have to be extended into the time dimension to account f o r unsteady flows along the r i v e r , adding a possibly unmanageable l e v e l of complexity to the solution method. The a p p l i c a t i o n of a dynamic programming formulation to t i d a l r i v e r water q u a l i t y problems poses an ad d i t i o n a l problem. A DP algorithm would be required which would optimize over upstream and downstream flows and BOD/DO int e r a c t i o n s (an o s c i l l a t i n g scheme). The development of such a technique would be a f r u i t f u l area of research. An ad d i t i o n a l order of complexity would ar i s e i n accounting f or the stochastic nature of the processes involved. 2.3.4 Stochastic Programming Models The a p p l i c a t i o n of LP, NLP and DP formulations to the water qua l i t y problem discussed previously were p r i m a r i l y d e t e r m i n i s t i c . The p r o b a b l i s t i c nature of the studies reviewed was not s p e c i f i c a l l y addressed. Random unit costs, resource constraints and unit resource usage could s i g n i f i c a n t l y complicate the optimizations schemes reviewed. Yeh (1985) documented the a p p l i c a t i o n of stochastic programming techniques i n the reso l u t i o n of rese r v o i r management problems. The a p p l i c a t i o n of these techniques to water q u a l i t y problems i s also growing. - 33 -Yeh (1985) discussed the a p p l i c a t i o n of chance-constrained LP (with and without recourse) and r e l i a b i l t y - c o n s t r a i n e d NLP and DP methods. Chance-constrained programming problems are of the form (as studied by Charnes, Cooper and Symonds, 1958) of: ^max/minj C_" X subject to Pr [ A X <_ B ] >_ e_ where C i s the cost vector 5C i s the decision vector A i s a deterministic matrix B i s a vector of random variables e i s a vector of exceedance p r o b a b i l i t i e s (Kolbin, 1976) Assuming normality and independence for the elements of 13, the above constraint can then be converted to a determininstic equivalent as: l-e) LJ " i — ~ i 1 f o r a l l i . Using t h i s type of formulation Lohani and Thanh (1978) addressed the minimization problem of t o t a l operating cost of BOD removal by optimizing the removal f r a c t i o n required by 'n' treatment plants to meet DO standards i n a common receiving r i v e r . Their problem was formulated n as: iTYUn i - t and <i'<r.<Ri' •, Y l ? o (2-17» - 34 -where i s the ultimate BOD waste discharge at reach i . B^ i s the allowable BOD to maintain DO i n reach i Y^ i s the unremoved f r a c t i o n of BOD i n reach i L i j """s t^ i e ^ e c o m P o s i t : L o n r a t i o for BOD between reaches i and j q^ i s the waste discharge flow at reach i i s the stream discharge at i a^ i s the r e l i a b i l i t y parameter where the constraints are the allowable r i s k of v i o l a t i n g the performance requirements (namely the DO standard) and the bounds on treatment f r a c t i o n with stream flow as the only random quantity. They used a steady-state model due to Camp i n e s t a b l i s h i n g the values of B. and L,.. By s e n s i t i v i t y analysis, these researchers obtained the minimum operating cost function for various r e l i a b i l i t y l e v e l s and various DO l i m i t s f o r the Hsintian River i n Taiwan. Jacoby and Loucks (1972) used a stochastic programming model to determine a small set of superior management plans for further evaluation with a simulation model i n t h e i r study of the Delaware River which focussed on the optimal mix of recreation, h y d r o e l e c t r i c and water supply p r o j e c t s . They also used a dynamic program to determine the optimal project completion dates. Burn and McBean (1985) formulated a program to minimize treatment costs at two e f f l u e n t sources on the Speed River near Guelph, Ontario which had six monitoring stations. They also assume steady-state conditions using the Streeter-Phelps model. Their problem was: max T" V B. W. d. . X. . 4 - i - 3 1 i ] i ] J c where B. i s a weighting factor - 35 -w. d.. i s the random variable which i s the combination of random pollutant loading and random transf e r c o e f f i c i e n t s of BOD at s i t e i monitorred at s i t e j X_ i s the removal e f f i c i e n c y (decision variable) of treatment at s i t e i on s i t e j subject to L e X * k a budgetary constraint where i s a piecewise l i n e a r cost curve and C L i s the budget. Minimum and maximum treatment l e v e l constraints were also used with a range of treatment 'equity'. The programming problem was solved as a dual LP and compared with r e s u l t s from a Monte Carlo simulation model with less than 5% variance i n the r e s u l t s . These researchers conclude that a programming model which accounts for the stochastic nature of variables involved would enhance the agreement between screening models and simulations allowing for better optimal solutions. 2.3.5 Control Methods The optimal a l l o c a t i o n of BOD can also be formulated as a c l a s s i c a l optimal control problem solvable by continuous or d i s c r e t e methods. These formulations address the problem of real-time control of BOD and DO i n a r i v e r by load regulation or a r t i f i c i a l aeration. Tarassov et a l (1969) used an algorithm akin to Pontryagin's maximum p r i n c i p l e based on the method of c h a r a c t e r i s t i c s . Using the Streeter-Phelps model as a "quasi-linear, f i r s t - o r d e r , m ultivariable, d i s t r i b u t e d parameter optimal control model" with pertubation techniques, they developed a Hamiltonian function and associated costate equations and boundary conditions. These equations were solved for a - 36 -quadratic objective function by numerical i n t e g r a t i o n generating various state response surfaces. They examined aerator controls i n time and space, space only and time only for an i n i t i a l l y p o l l u t e d r i v e r . Davidson and Bradshaw (1970) studied the same problem using the maximum p r i n c i p l e and noted that the solution of the problem was also possible using dynamic programming techniques. However, t h e i r objective was to develop a mathematically t r a c t a b l e model for numerical s o l u t i o n which tested the s e n s i t i v i t y of the optimal control p o l i c y to stream temperature. Tabak and Kuo (1970) presented extensive discussion on the re l a t i o n s h i p s between mathematical programming methods and solution methods for optimal control problems showing t h e i r complementary and interchangeable uses. These two early works were assumed steady-state BOD/DO processes while l a t e r studies surveyed incorporated dynamic BOD/DO models with steady hydrodynamics. Tamura (1974) used the Streeter-Phelps model to develop multi-dimensional, high-order difference equations which accounted for BOD/DO transport delays. This formulation was solved numerically using d u a l i t y and decomposition techniques i n minimizing the deviation between the system output DO, due to BOD loadings, and regulatory DO concentration l e v e l s . Gourishankar and Lawson (1978) used a steady-state model to study a multivariable control formulation of the BOD p o l l u t i o n problem. They developed control equations to tes t a BOD dumping c o n t r o l l e r and one f o r both dumping and instream aeration. Ozgaren et a l (1975) developed a control strategy f o r aeration i n a r i v e r system where observations are time-lagged while Olgac et a l (1976) studied the optimal a l l o c a t i o n of measurement and control resources by introducing white noise processes into the BOD/DO - 37 -demand and DO measurements. Ja i n and Denn (1976) applied control techniques to short-term BOD impulses i n the Delaware Estuary and concluded that t h e i r optimal p o l i c y required rapid control responses that were impossible to achieve i n p r a c t i s e . Gourishankar and Ramar (1977) designed a continuous DO c o n t r o l l e r f or a dynamic BOD/DO non-tidal steady discharge r i v e r model which was enhanced by Kudva and Gourishankar (1977) who considered the case of an observer of the same system where inputs were i n a c c e s s i b l e . Gourishankar and Lawal (1978) presented a d i g i t a l DO c o n t r o l l e r design with lower instrumentation costs for equivalent performance compared to previously designed continuous c o n t r o l l e r s . Spear and Hornberger (1983) discussed the impacts of parameter uncertainty on c o n t r o l l e r operation showing that the simple c o n t r o l l e r studied could control DO l e v e l s eighty-four percent of the time for the River Cam, England. Additional to these works on real-time c o n t r o l l e r design, were studies applying m u l t i l e v e l of h i e r a r c h i c a l techniques to real-time and seasonal control s t r a t e g i e s . For example, Haimes (1976) presented an a n a l y t i c two-level formulation of a regional p o l l u t i o n problem using the Streeter-Phelps model. The f i r s t l e v e l represented l o c a l p o l l u t e r s on the r i v e r with t h e i r respective treatment costs while the second l e v e l represented a r i v e r basin wide regional authority with respective s o c i a l costs. The r e s u l t i n g NLP problem was solved using gradient p r o j e c t i o n methods. A t h r e e - l e v e l scheme was also proposed where the regional a u t h o r i t i e s would report to a central one. Singh (1977) used h i e r a r c h i c a l decomposition techniques i n the optimal control of p o l l u t a n t loading so as to return the state of the system to 'steady' conditions at minimum cost, also for the River Cam. Chichester (1976) - 38 -also presented a m u l t i - l e v e l formulation based on a compartment dynamic p o l l u t a n t model of B e l l a and Dobbins (1968). Three control methods were discussed: aeration, dumping and a combination of the two. The control problem was s p l i t into sub-problems based on the s o l u t i o n of the state, control and costate equations derived for a one-dimensional t i d a l r i v e r which were combined h i e r a r c h i c a l l y to a 'co-ordinator'. The numerical solution for the optimal real-time control of DO required computing resources of 9 MB of core and 45 seconds CPU time for a four reach t i d a l r i v e r with s p a c i a l l y c o n t r o l l e d aeration and p o l l u t a n t dumping. Thus the a p p l i c a t i o n of m u l t i - l e v e l techniques appear r e s t r i c t e d by the number of sub-problems that can be managed computationally and conceptually. 2.4 Summary The management of r i v e r systems i s d i f f i c u l t and complex process aided by the judicious choice and use of well developed and v a l i d a t e d models. Friedman et a l (1984) i n an i n v e s t i g a t i o n of the p o t e n t i a l uses of models i n the U.S., states that the use "of mathematical models have s i g n i f i c a n t l y expanded the nation's a b i l i t y to understand and manage i t s water resources" and that "models have the p o t e n t i a l to provide even greater benefits for water resource decision making i n the future". The stochastic nature of the water q u a l i t y processes also need to be e x p l i c i t l y addressed by these models (Ward and L o f t i s , 1983) while the combined use with optimization models would guide managers and users towards more optimal resource use (Loucks, 1983). 3 . MODEL DEVELOPMENT T h i s c h a p t e r p r e s e n t s t h e f o r m u l a t i o n o f BOD c o n t r o l d i s c u s s e d i n c h a p t e r 1 as a m a t h e m a t i c a l programming p r o b l e m , f o l l o w e d by the d e v e l o p m e n t o f a w a t e r q u a l i t y m o d e l u s e d i n t h i s s t u d y . S p e c i f i c a l l y , a d y n a m i c BOD/DO model and a d i r e c t s e a r c h a l g o r i t h m a r e d e s c r i b e d . 3 . 1 P r o b l e m F o r m u l a t i o n T h e p r o b l e m o u t l i n e d i n c h a p t e r 1 i s one o f m a x i m i z i n g BOD l o a d i n g s o f a r i v e r s y s t e m w h i l e m a i n t a i n i n g a' c e r t a i n l e v e l o f DO f o r f i s h e r y r e q u i r e m e n t s . S i n c e t h e r i v e r s y s t e m i s b o t h d y n a m i c and s t o c h a s t i c , an a b s o l u t e DO s t a n d a r d w o u l d be d i f f i c u l t t o impose due t o i t s r i g i d i t y as compared t o one w h i c h a l l o w s an a c c e p t a b l e f r e q u e n c y o f n o n - c o m p l i a n c e . The c o n c e p t o f f l e x i b l e w a t e r q u a l i t y s t a n d a r d s has been u s e d more r e c e n t l y f o r m u l t i - u s e r w a t e r s h e d s i n t h e U n i t e d S t a t e s and the E u r o p e a n Economic Community ( B e c k , 1985 ) . 3 . 1 . 1 O b j e c t i v e A v a r i e t y o f o b j e c t i v e f u n c t i o n f o r m s a r e p o s s i b l e f o r u s e a s s u r v e y e d i n s e c t i o n 2 . 3 . The o b j e c t i v e as s t a t e d above c o u l d be c o n s t r u c t e d a s : a) a u t i l i t y f u n c t i o n b a s e d o n , s a y , s o c i a l c o s t s o r p e r h a p s agency and u s e r p r e f e r e n c e c u r v e s ; o r b ) a c u m u l a t i v e f u n c t i o n o f BOD l o a d i n g s . T h e s e f u n c t i o n s c o u l d be e i t h e r l i n e a r o r n o n - l i n e a r , c o n v e x or c o n c a v e . T h e f o r m b) above was c h o s e n p r i m a r i l y f o r i t s s i m p l i c i t y and a m e n a b i l i t y t o l i n e a r f o r m u l a t i o n s . Form a ) , once d e f i n e d , c o u l d a l s o be - 40 -e a s i l y f o r m u l a t e d . A s i m i l a r c h o i c e a l l o w e d G o r e l i c k (1982) t o f o r m u l a t e a s i m p l e l i n e a r program f o r a g r o u n d w a t e r management p r o b l e m . BOD l o a d i n g s were assumed t o o c c u r o v e r d i s c r e t i z e d r i v e r r e a c h e s w h i c h were j u d g e d t o be s u f f i c i e n t l y f i n e so as t o g e n e r a t e m e a n i n g f u l management p o l i c i e s . T h e s e r e a c h e s were assumed t o c o n t a i n b o t h p o l l u t a n t l o a d i n g s and m o n i t o r i n g s t a t i o n s . The i s s u e o f t h e o p t i m a l l e v e l o f r i v e r d i s c r e t i z a t i o n was n o t a d d r e s s e d . T h u s t h e o b j e c t i v e f u n c t i o n was s i m p l y max x i w i ( 3 . 1 ) U i w h e r e : x^ i s t h e c o n t i n u o u s , s t e a d y n o n - p o i n t BOD l o a d i n g a l o n g r e a c h i o f t h e r i v e r ; [M.L - - ' - ] w± i s t h e l e n g t h o f r e a c h i ; [L] N i s t h e number o f r e a c h e s i n t h e r i v e r . (Note t h a t w^ c o u l d be c h o s e n t o r e f l e c t r e l a t i v e i m p o r t a n c e or e c o n o m i c v a l u e o f t h e d i f f e r e n t r i v e r r e a c h e s . ) 3 . 1 . 2 C o n s t r a i n t s T h e c o n s t r a i n t on t h e m a x i m i z a t i o n p r o b l e m was t h a t o f a d h e r a n c e t o a r e g u l a t o r y DO s t a n d a r d f o r a l l r e a c h e s o f t h e r i v e r o v e r a d e c i s i o n h o r i z o n , w i t h some o v e r a l l l e v e l o f c o m p l i a n c e . T h u s i n d i v i d u a l r e a c h r e g u l a t i o n s were n o t assumed n o r was the o v e r a l l r i s k a p p o r t i o n e d by r e a c h . T h i s may have been p o s s i b l e w i t h a p r i o r i i n f o r m a t i o n on s e n s i t i v e s i t e s a l o n g t h e r i v e r where s p e c i f i c f l e x i b l e s t a n d a r d s may be u s e f u l f o r p o l l u t a n t c o n t r o l . The r e g u l a t o r y s t a n d a r d used was t h e minimum i n - s t r e a m DO l e v e l s f o r r e s i d e n t and m i g r a t o r y f i s h s p e c i e s , w h i c h a r e p r i m a r i l y s a l m o n i d s f o r P a c i f i c t i d a l r i v e r s . - 41 -The c o n s t r a i n t s were f o r m u l a t e d as n o n l i n e a r e q u a l i t y and i n e q u a l i t y e q u a t i o n s : ( ! ) y i j t = * j t ( x i + x o i » I D B > (2) P j t = _ > > i j t V (3) z A j = g j t C P j t , 1 1 ^ ) ) (4) u j t = r 1 i f z j t > C r e g [ 0 i f z j t < C r e g (5) S = ( 1 / N t ) £ £ u j t > ^ > 0 (6) x ± > 0 } fc Y i j t > o z j t > o where i , j = 1 , . . , N r e a c h e s , i b e i n g t h e i n p u t r e a c h a n d j t h e o u t p u t r e a c h and t = 1 , . . , X , t h e d e c i s i o n h o r i z o n w i t h A t d e t e r m i n e d f r o m s t a b i l i t y c r i t e r i a o f t h e f i n i t e - d i f f e r e n c e m o d e l u s e d ( s e c t i o n 3 . 2 ) . O f t h e c o n s t r a i n t s a b o v e : (1) a r e t h e d i f f e r e n c e e q u a t i o n s o f t h e BOD m o d e l ( s e c t i o n 3 . 2 . 2 . 3 ) ; (2) r e p r e s e n t the a d d i t i v i t y a s s u m p t i o n s whereby t h e BOD c o n c e n t r a t i o n s o f s i m u l t a n e o u s l o a d i n g s i n m u l t i p l e r e a c h e s a r e e q u i v a l e n t t o t h e sum o f t h e e f f e c t s o f i n d i v i d u a l l o a d i n g s ; (3) a r e t h e d i f f e r e n c e e q u a t i o n s o f t h e DO m o d e l s ( s e c t i o n 3 . 2 . 2 . 4 ) ; (4) and (5) c o m p r i s e t h e p r o p o r t i o n a l i t y c o n s t r a i n t . I n t h e above c o n s t r a i n t s : y i j t * s t n e B ^ ) u c o n c e n t r a t i o n i n r e a c h j a t t i m e t due t o l o a d i n g i n r e a c h i a t t i m e t - 1 ; P j t i s t h e t o t a l BOD c o n c e n t r a t i o n i n r e a c h j a t t i m e t due t o t h e t o t a l l o a d i n g s x^ i n r e a c h e s i • = ! , . . , N a t t i m e t - 1 ; - 42 -Z j t i s the DO c o n c e n t r a t i o n i n r e a c h j a t t i m e t due t o the BOD c o n c e n t r a t i o n i n r e a c h j ; U j t i s an i n d i c a t o r v a r i a b l e o f Z j t v e r s u s C r e g f o r r e a c h j a t t i m e t ; S i s t h e p r o p o r t i o n o f U j t f o r j = l , . . , N and t = l , . . , t w h i c h meet the C r e g r e g u l a t i o n l i m i t ; mg, mj) a r e m a t r i c e s o f m o d e l p a r a m e t e r s ; ^ i j t ( x i » m B ) * s t h e s e t o f e q u a t i o n s 3 . 1 1 - 3 . 1 3 , 3 . 1 7 - 3 . 1 8 ; 8 j t ( P j t » m D ) * s t h e s e t o f e q u a t i o n s 3 . 1 1 - 3 . 1 3 a n d 3 . 1 9 ; X Q ^ i s t h e i n i t i a l BOD c o n c e n t r a t i o n s i n r e a c h i a t t=0; C r e g i s t h e DO l i m i t f o r a l l r e a c h e s ; i s t h e a c c e p t a b l e c o m p l i a n c e l e v e l ; a n d N i s the number o f r e a c h e s . C o n s t r a i n t s (4) a n d (5) c a n be v i e w e d i n t h e s p i r i t o f t h e c l a s s i c a l j o i n t - c h a n c e c o n s t r a i n t w h i c h d e s c r i b e t h e p r o b a b i l i t i e s o f e x c e e d a n c e o f s t o c h a s t i c v a r i a b l e s . T h e p r o b l e m s t u d i e d h e r e does n o t a d d r e s s t h e s t o c h a s t i c n a t u r e o f the p r e v i o u s l y d e f i n e d v a r i a b l e s b u t d e a l s w i t h the a v e r a g e s t a t i s t i c s o f a d e t e r m i n i s t i c p r o c e s s . Had t h e s t o c h a s t i c n a t u r e o f the p o l l u t a n t and h y d r o d y n a m i c p r o c e s s e s mode l e d been i n c l u d e d i n the n o n - l i n e a r p r o g r a m , a n a l y t i c f a c i l i t a t i o n o f t h e s o l u t i o n w o u l d have b e e n h i g h l y u n l i k e l y . Wagner and M i l l e r (1965) d i s c u s s e d t h e c l a s s i c a l j o i n t - c h a n c e c o n s t r a i n t p r o b l e m w i t h random l i n e a r c o e f f i c i e n t s : max z = £ X s u c h t h a t A X < b and X >_ £ ( 3 . 3 ) and I T P r [ M a ; X ] >, ^ ^ O , < 3 - 4 ) S b e i n g t h e s e t o f chance c o n s t r a i n t s ) where a^ i s a random v e c t o r and b± a r e known. T h e y , h o w e v e r , c o n s i d e r e d a f o r e a c h i a s m u l t i n o r m a l , w i t h mean a n d c o v a r i a n c e m a t r i x v , and as c o m p l e t e l y i n d e p e n d e n t . - 43 -A n a l y t i c a l s o l u t i o n t e c h n i q u e s p r o p o s e d by P r e k o p a ( S e n g u p t a , 1972) and S e n g u p t a (1972) a l s o r e q u i r e ^ n i c e ^ d i s t r i b u t i o n a l a s s u m p t i o n s . T h e s e r e s t r i c t i o n s a r e t o o n a r r o w f o r t h e s t o c h a s t i c v e r s i o n o f t h e p r o b l e m a t hand where a± a r e n o n - n o r m a l w i t h s p a c i a l a n d t e m p o r a l a u t o -c o r r e l a t i o n . T h e s i m p l i f i e d d e t e r m i n i s t i c a p p r o x i m a t i o n method o f c o n s t r a i n t s (4) and (5) a b o v e , b a s e d on a v e r a g e d s t a t i s t i c s f o r j = l , . . , N and t=l , . . ,"C may n o t be an e l e g a n t s o l u t i o n m e t h o d , w i t h p o s s i b l e a n a l y t i c a l s i m p l i f i c a t i o n s , but a p p e a r e d v i a b l e w i t h t h e n u m e r i c a l t e c h n i q u e s and p o l l u t a n t m o d e l s p r e s e n t e d l a t e r i n t h i s c h a p t e r . T h e r e i s one i m p o r t a n t d i f f e r e n c e between f u n c t i o n s g j t and f i j t w h i c h a l l o w s one t o s i m p l i f y t h e n u m e r i c a l method o f s o l u t i o n o f t h e o p t i m i z a t i o n p r o b l e m as f o r m u l a t e d . T h e BOD f u n c t i o n f i j t I s b a s e d on a l i n e a r lumped p a r a m e t e r s o l u t e t r a n s p o r t m o d e l . A l i n e a r i n c r e a s e i n i n p u t l o a d i n g s r e s u l t s i n an e q u i v a l e n t l i n e a r i n c r e a s e i n o u t p u t c o n c e n t r a t i o n s . T h i s p r o p e r t y was e x p l o i t e d by G o r e l i c k (1982) i n d e v e l o p i n g a s i m u l a t i o n - b a s e d l i n e a r p r o g r a m i n m a n a g i n g g r o u n d w a t e r p o l l u t i o n s o u r c e s . T h i s w o u l d t h e a l l o w f i j t t o be r e p l a c e d by a l i n e a r f u n c t i o n o f u n i t l o a d i n g s o r fijtCxi+Xo 1^) = ( x i / u i ) * f i j t ( u i » m B ) + x o i where u ^ a r e u n i t l o a d i n g s i n r e a c h i . T h i s a l s o a l l o w e d f o r the a d d i t i v e a s s u m p t i o n o f c o n s t r a i n t ( 2 ) . T h e f u n c t i o n g j t i s a n o n - l i n e a r DO model and t h u s c a n n o t be c o n v e n i e n t l y s i m p l i f i e d . T h u s , w i t h i n t h e s o l u t i o n p r o c e d u r e ( s e c t i o n 3 . 3 . 1 ) t h e BOD a n d DO f u n c t i o n s were s e p a r a t e d , t h e l a t t e r b e i n g d i r e c t l y i n c o r p o r a t e d i n t o t h e o p t i m i z a t i o n a l g o r i t h m o f s e c t i o n 3 . 3 . 3 . T h e s p e c i f i c a t i o n o f t h e a c c e p t a b l e l e v e l o f c o m p l i a n c e i s a d i f f i c u l t p r o c e s s and may n e e d t o be r a t i o n a l i z e d u s i n g f i s h p o p u l a t i o n - 44 -m o d e l s . The c o m p l i a n c e l e v e l , r i v e r d i s c h a r g e , r e g u l a t o r y D O , BOD d e c a y c o e f f i c i e n t and l o n g i t u d i n a l d i s p e r s i o n were u s e d p a r a m e t r i c a l l y i n p o l i c y s e n s i t i v i t y a n a l y s i s . 3 . 2 W a t e r Q u a l i t y M o d e l T h e w a t e r q u a l i t y m o d e l u s e d i n t h i s s t u d y c o m p r i s e d two c o m p o n e n t s : I ) a h y d r o d y n a m i c r i v e r m o d e l ; I I ) an u n s t e a d y BOD/DO m o d e l . M o d e l I s i m u l a t e d r i v e r s t a g e and d i s c h a r g e d y n a m i c s w h i l e m o d e l I I p r e d i c t e d DO and BOD c o n c e n t r a t i o n s due t o BOD l o a d i n g s . The u n c o u p l i n g o f t h e w a t e r q u a l i t y m o d e l i n t o t h e s e two components was b a s e d on t h e a s s u m p t i o n t h a t t h e two phenomena were w e a k l y i n t e r d e p e n d e n t s u c h t h a t t h e d e n s i t y o f t h e r i v e r w a t e r s was n o t a f f e c t e d by t e m p e r a t u r e , s a l i n i t y and h y d r o - b i o l o g i c a l e f f e c t s ( G r o m i e c e t a l , 1 9 8 3 ) . T h i s t h e n a l l o w e d : a) s i m p l i f i c a t i o n o f model f o r m u l a t i o n and d e v e l o p m e n t ; b) t h e i n c o r p o r a t i o n o f a h y d r o d y n a m i c model u s e d i n a p r e v i o u s m o d e l i n g s t u d y ( J a m a l , 1980) ; c ) the d e v e l o p m e n t o f a s e p a r a t e BOD/DO m o d e l ; d) i n d e p e n d e n t v a l i d a t i o n o f m o d e l s I a n d I I ( w h i c h t h o u g h n o t as r o b u s t a p r o c e d u r e as combined v a l i d a t i o n , was n e c e s s a r y due t o t h e l a c k o f d a t a ) . 3 . 2 . 1 M o d e l I : H y d r o d y n a m i c R i v e r M o d e l 3 . 2 . 1 . 1 B a s i c A s s u m p t i o n s T h e m o d e l was d e v e l o p e d f o r m o d e l i n g o f t h e r i v e r s y s t e m and d i d n o t i n c l u d e p r o c e s s e s r e l a t e d t o g r o u n d w a t e r r e c h a r g e and d i s c h a r g e , - 45 -overland flow, evapotranspiration, small t r i b u t a r y inflows, s a l i n i t y e f f e c t s , wind e f f e c t s , bed load and sediment e f f e c t s . The impact of licensed i r r i g a t i o n uptakes was reviewed and i s discussed i n section 4.4. It was also assumed that the hydrodynamics of the r i v e r could be adequately described by the set of one-dimensional d i f f e r e n t i a l equations of continuity and momentum. 3.2.1.2 Basic Equations Continuity: (3.6) ^5 ^ vs/Db - o Ox r > t 3 -1 where Q i s the r i v e r discharge [L .T ] W i s the bank width [L] h i s the r i v e r stage to (Geologic Survey of Canada) datum [L] x i s the distance along the r i v e r axis [L] t i s the time dimension [T] Momentum: CHJ ^ = d h - 9 M O (3.7) Ofc Ox £>x cM where u i s the r i v e r v e l o c i t y [L.T _2 g i s the acceleration due to gravity [L.T ] c i s the Chezy c o e f f i c i e n t d i s the depth of flow [L] 2 A i s the cross-sectional area [L ] (Gromiec et al,1983) Substitution of the averaged discharge as Q=vA i n equation 3.7 and with i n - s u b s t i t u t i o n of equation 3.6 leads to A O t A*5W Alttik A l 5 x i ~C?&& (3'8) - 46 -The boundary conditions for equations 3.6 and 3.8 were the upstream inflow discharges into the r i v e r system and downstream t i d a l heights. I n i t i a l conditions were the stage and discharge at various cross-sections along the r i v e r . In t h i s study two f i n i t e - d i f f e r e n c e schemes were attempted i n the so l u t i o n of equations 3.6 and 3.8 subject to boundary conditions. The f i r s t attempt was the development of an i m p l i c i t difference scheme based on Dronkers (1969). This scheme required the so l u t i o n of a set of simultaneous equations for cross-sections along the r i v e r at each time step. The 'double sweep' method for matrix so l u t i o n was used with recorded i n i t i a l and boundary conditions. The method, however, f a i l e d to adequately simulate the dynamics of a t i d a l gate which existed on the r i v e r of concern and had a primary influence on the r i v e r flows. The method was abandoned i n favour of an e x p l i c i t scheme presented by Ages (1973). An e x p l i c i t f i n i t e difference formulation of equations 3.6 and 3.8 was presented by Jamal (1980) using a 'leap frog' scheme and was adapted to f i t the BOD/DO model discussed i n section 3.2.3. This formulation i s summarized below. The r i v e r was schematized as shown i n figure 3.1 with the x - d i r e c t i o n p o s i t i v e upstream. Equations 3.6 and 3.8 were d i s c r e t i z e d by c e n t r a l differences to y i e l d equations 3.9 and 3.10 i n a scheme represented i n f i g u r e 3.2 showing i n i t i a l and boundary conditions f o r the Nicomekl River i n Surrey, B.C. Details of the development of equations 3.9 and 3.10 are found i n Jamal (1980). The t i d a l gate was simulated by simply f o r c i n g to zero the discharge at the gate for a l l periods that the r i v e r discharge was i n the p o s i t i v e upstream d i r e c t i o n . The c a l i b r a t i o n parameter was the Manning's 'n' number which r e f l e c t s - 48 -m m+ 1 < + 2 < Ax + Ax _Ll m m+1 (3.9) ,n+l m R L A " + 1 + A V m n + 1 m+ 1 A x +Ax , m m-1 T < 1 •CI> A n + l + A n + l _ A n _ An m+1 m-1 m+1 m-1 At / . n + l . A n + l x 2 (Am + A m + 1 ) w n + l + w n + l H n + 1 + H n + 1 fln n m m + 1 . m+1 m-1 m-1 m+1 .»< <Ci-C;> / A n + l . . n + l v 2 2 , u n + l „ 0 . (A +A .) c (H -H ) m m+1 m m m L 1 1 | 1 • • Boundary Conditions (Q below Murray Creek at 206th S t * Langley) - 6V -- 50 -the roughness of the r i v e r channel. The numerical s t a b i l i t y of the scheme i s based on a 'At' chosen from the well known c r i t e r i a An overview flowchart"of t h i s model i s presented i n figu r e 3.3 where the stage and discharges are updated at each time step for a l l r i v e r cross-sections. 3.2.3 Model I I ; The BOD/DO model The models surveyed i n section 2.2.2 were of three types: a) models based on PDE's and solved by f i n i t e difference equations b) as a) above but solved by c h a r a c t e r i s t i c methods; c) compartment or box models representing mass balances on r i v e r reaches. The development of t h i s model was based of the study of B e l l a and Dobbins (1968) (type c) above) which, though not state-of-the-art (Ambrose, 1986), was of a simpler formulation than a) or b) and thus easier to incorporate into the numerical hydrodynamic model of section 3.2.2. 3.2.3.1 Basic Assumptions The water q u a l i t y model developed was based s t r i c t l y on r i v e r BOD/DO dynamics. No attempt was made to develop e c o l o g i c a l model components due to the complexity of such a task e s p e c i a l l y with such paucity of data and lack of f i e l d resources. The model also does not simulate overland BOD loadings due to storm events or f a l l runoff, such as Tubbs and Haith (1981), although these have been recognized as s i g n i f i c a n t p o l l u t a n t processes (Moore, 1985). It i s assumed that the BOD decay was f i r s t - s t a g e carbonaceous process with no second-stage nitrogenous influence. This was a si m p l i f y i n g assumption since - 51 -R E A D INITIAL A N D B O U N D A R Y C O N D I T I O N S R E A D C H A N N E L D A T A INCREMENT TIME BY Ate-INCREMENT ODD NUMBERED CROSS-SECTIONS c-FROM EQUATION OF CONTINUITY COMPUTE HEIGHT AT SECTION m Z Z J L COMPUTE AVERAGE CHANNEL PARAMETERS INCREMENT EVEN NUMBERED CROSS-SECTIONS C~~A FROM EQUATION OF MOTION COMPUTE DISCHARGE AT SECTION m END OF TIME PERIOD ? NO OUTPUT OF SIMULATED HEIGHTS AND DISCHARGES Figure 3.3: Flow Chart of the Hydrodynamic Model (Jamal, 1980) - 52 -second-stage processes are often present (Ambrose, 1986). It was also assumed that oxygen sources and sinks i n the form of d e t r i t u s , benthos and algae could be accounted for by some base l e v e l consistent with f i e l d experience. Diurnal photosynthetic models are a v a i l a b l e i n the l i t e r a t u r e (e.g. Jorgensen, 1983) but were not used due to the lack of data for v a l i d a t i o n . Steady net oxygen rates were used to account for photosysthesis and r e s p i r a t i o n . The stream temperature was also assumed constant for the time period of i n t e r e s t such that temperature dependent parameters were non-varying. It was assumed that the processes of convection, dispersion and decay (or oxidation) were s u f f i c i e n t to describe the one-dimensional BOD/DO dynamics with homogeneous concentrations i n each reach. 3.2.3.2 Basic Equations Convection: the transport of the pollutant by the movement of the r i v e r waters downstream or upstream. The mass balance of a pollutant over a stream segment for a time i n t e r v a l i s described as: mass at end mass at s t a r t mass convected mass convected of i n t e r v a l = of i n t e r v a l + i n during - out during i n t e r v a l i n t e r v a l or (for downstream flows), L(n,T+4T)A(n,T+AT)Ax n = L(n,t)A(n,T)Ax n + UA(n+1/2,T)L(n+1,T)AT - UA(n-1/2,T)L(n,T)AT (3.11) which can be simply solved for L(n,T+AT) where: L i s the pollutant concentration [M.L - 3] A i s the r i v e r c r oss-sectional area Ax, n i s the length of reach n (+ve upstream) where n i s located at even numbered sections - 53 -of f i g u r e 3.2. [ L ] AT i s the time i n t e r v a l for T=1, . . , T [ T ] U i s the r i v e r v e l o c i t y (+ve upstream) [ L . T ~ 1 ] 3 — 1 UA i s the r i v e r discharge (+ve upstream) [ L , T ] n i s the midpoint of a reach A * n running from n-1/2 to n+1/2 as i n fi g u r e 3.4, n=2,4,..,N. For upstream flows the mass balance equation i s : L(n , T+AT)A(n , T+AT)Ax n = L ( n , T ) A ( n , T ) A x n + UA(n-1 / 2 , T ) L(n - 1 , T)AT - UA(n+1 /2 , T ) L(n , T)AT (3.12) B e l l a and Dobbins (1968) discussed to some length the d e f i n i t i o n s of mixing and dispersion to ensure that the convection equations above c o r r e c t l y represent the physical processes involved. Dispersion: the d i f f u s i o n of the pollutant due the mixing and exchange of r i v e r waters while the r i v e r i s flowing. The mass balance on a stream segment i s : L (n ,T+AT)A(n,T) = DA(n-1/2 ,T)AT [ L ( n - 1 , T ) - L ( n , T ) ] + L ( n , T ) A ( n , T ) + DA(n+1 /2 ,T)AT [L(n+1 ,T) - L ( n , T ) ] A x 2 n (3.13) where D i s the dispersion c o e f f i c i e n t for the r i v e r as a whole or by reach. In f i n i t e - d i f f e r e n c e modeling, pseudo- (or numerical) dispersion often r e s u l t s due to the numerical approximations made. This dispersion can be minimized, as suggested by B e l l a and Dobbins (1968), by reducing the c o e f f i c i e n t D by a corresponding amount or by choosing Ax,^ and AT appropriately. The pseudo-dispersion was estimated as D A(n+1/2 ,T) = UA(n+1/2 ,T)[Ax - UA(n+1/2,T)AT ] (3.14) P 2 A(n , T + A T ) such that D i s corrected as D' = D - D . A l t e r n a t i v e l y , Ax and A T - 54 -DOWNSTREAM T + A T T + 2 A T Figure 3.4: Convection of a Slug Load (Bella and Dobbins, 1968) - 55 -can be chosen from ( UA(n+1/2,T) ) = Ax R (3.15) ( A(n,T+AT) )max AT I f , however, i s 'small' compared to D then i t can be neglected. In t h i s study the f i r s t method was used. The s t a b i l i t y c r i t e r i a for the dispersion model i s given as ( D A T ) < 1 (3.16) ( A x n 2 ) 2 to prevent o s c i l l a t i n g errors. Decay: i s the process whereby the pollutant i s biochemically neutralized within the r i v e r . This was represented i n a f i r s t - o r d e r decay equation as: L(n,T+AT) a L(n,T) - K.,(n,T)AT [(1-0)L(n,T) + ©L(n,T+AT)] (3.17) where K^(n,T) i s the decay c o e f f i c i e n t which may vary i n time and space, though i s was considered constant i n t h i s study. © i s a c o e f f i c i e n t between 0 and 1 which averages the pollutant over a simulated time i n t e r v a l , AT. A value of 0.5 was chosen f o r simple averaging. Addition: pollutant sources can be added either as impulse or continuous loads over a reach of the r i v e r . A continuous load i s represented as a serie s of impulse loadings over each time i n t e r v a l . This was modeled as L(n,T+AT) = L(n,T) + m(n,T)AT (3.18) A(n,T)Ax n where m(n,T) i s the average rate of pollutant added to reach n over AT and can vary i n any reasonable manner. 3.2.3.3 BOD model: Equations 3.<J| to 3.J& above are s u f f i c i e n t to model the fate of BOD loadings i n a t i d a l r i v e r and as such could be c a l i b r a t e d and - 56 -v a l i d a t e d separately from the DO model discussed i n section 5.2. The model was used i n a "multi-step" scheme proposed by B e l l a and Dobbins (1968) where for every time step, and a l l reaches, the pollutant's f a t e i s simulated by sequentially modeling the processes of addition, convection, dipersion and decay i n an "add-con-disp-dk" procedure. This scheme i s flowcharted i n figure 3.5. 3.2.3.4 DO model Equations 3.11 to 3.16 and 3.18 can also be used to describe the processes involved i n the addition, convection and dispersion of oxygen. However, a source/sink DO model based on the mass balance concepts used e a r l i e r was required and presented as: C(n,T+AT) = C(n,T) - k 1(n,T)[(1-e)L(n,T)-0L(n,T+AT)]AT + k 2(n,T)[(1-6)(C s(n,T)-C(n,T)) + 6(C(n,T+AT)-C(n,T+AT) )]AT + 0(n,T)AT (3.19) where C(n,T) i s the DO concentration —3 i n reach n at time T [M.L ] k 2(n,T) i s the c o e f f i c i e n t of aeration [T~ 1] _3 C g(n,T) i s the saturated DO i n time and space [M.L ] — 1 0(n,T) i s the net of oxygen sources and sinks [M.T ] and the second term on the RHS of equation 3.22 was the e f f e c t of BOD concentration on DO l e v e l s . It i s t h i s term which makes the o v e r a l l DO model a non-linear lumped parameter model. For t h i s study, k 2 was computed from a reaeration equation due to Dobbins and 0'Conner (Covar, 1976) while the saturated DO l e v e l was taken as a constant based on f i e l d data and 0(n,T) was set at some minimum l e v e l due to benthic demand and photosynthesis. - 57 -INITIAL AND BOUNDARY CONDITIONS D, ki, k2, U, A, A X P , AT, —I z z : ADD: LaddCn,!) = UnJ) + ADDed BOD [Eqn 3.17] CON: Lcon(nJ) = CON{ Ladd(n/0 } [Eqn 3.11, 3.12] DISP: Ldisp(n.T) - DISP{ Lcon(nft) } [Eqn 3.13] DKY: Ldky(n,"T) = DKY{ Ldisp(nJ) } [Eqn 3.16] L(n,T+AT) o Ldky(n,"0 INCREMENT BY *T Figure 3.5: Flow Chart of the BOD Multi-Step Procedure - 58 -The combined BOD/DO model was used i n the multistep procedure i n f i g u r e 3.5 where each time step i t e r a t i o n was for a l l reaches i n the r i v e r . The case of a r t i f i c i a l r i v e r aeration could e a s i l y modeled by "adding' oxygen i n the multistep procedure and would be an i n t e r e s t i n g area for further research and was not addressed i n t h i s study though i t i s c urrently of great i n t e r e s t (Mavinic, 1986). 3.2.4 Combination of Water Quantity and Quality Models Since the interdependence of these two models was assumed to be weak the models developed previously can be coupled as a t y p i c a l linked water q u a l i t y model (Beck ,1985). Thus, the predictions of c r o s s - s e c t i o n a l area and r i v e r v e l o c i t y made by using the hydrodynamic model would be inputs to the water q u a l i t y model for the p r e d i c t i o n of BOD and DO with no feedback. This procedure i s flowcharted i n figure 3.6. 3.2.5 I n i t i a l and Boundary Conditions The i n i t i a l BOD and DO concentrations were set at p r i s t i n e r i v e r conditions of zero and one-hundred percent saturation r e s p e c t i v e l y i n the scheme of figure 3.7. These conditions would l i k e l y bias the predictions i f a suitable s t a b i l i z i n g period had not been incorporated into the simulation runs and a t y p i c a l system state not reached, e s p e c i a l l y for low flow conditions where the t i d a l gate may remain closed for extended periods of time. The upstream boundary conditions were also set as p r i s t i n e with no BOD entering the system at the uptream end while DO entered at saturated l e v e l s . Thus the upstream BOD loadings are assumed to have n e g l i g i b l e impacts on the r i v e r reaches of i n t e r e s t . The downstream conditions were more d i f f i c u l t to formulate - 59 -INITIAL BOUNDARY i PARAMETER VALUES. SIMULATE STAGE, DI5CHARGE. USING HYDRODYNAMIC MODEL FOR A t ("fig 3-3) I N C R E M E N T BY A t AT FROM. LAST CALL TO 50P/DO MODEL? TIME HORIZON ENDED ? STOP SIMULATE CONCENTRATIONS USING BOD/DO MODBL FOR &.J (Fig 35)N Figure 3.6: Flowchart of the Coupled Hydrodynamic and BOD/DO Models Figure 3.7: BOD/DO Model Schematic of River Reaches B O D / D O M O D E L PREDICTIONS REACH 1 REACU 2 REACH 3 REACH 4 s+e >)< >K 9H DISTANCE , Oft) 0 800 1600 7400 13200 30800 30600 3M00 40bOO 43800 4705O 5)010 UPSTREAM , I I I I I I I I I I I r ^ P P O T ^ T H A M I C M O D E L P R E D I C T I O N S "• 6) Q Q Q CRce&-secTION: 1 2 1 i TlXAL TIWO. KXJNDW GATE lO 11 12 1 WSG STATION Q = DBOtARGE M » HT- OF FLOW BELOW MURRAY CREEK - 61 -due to the t i d a l conditions. Two possible assumptions were available i n the absence of e s t u a r i a l q u a l i t y data: a) r i v e r and t i d a l d i l u t i o n i s s u f f i c i e n t l y great to d i l u t e the BOD and DO to n e g l i g i b l e l e v e l s ; b) the d i l u t i o n i s not s u f f i c i e n t to assume n e g l i g i b l e l e v e l s ; The former was assumed for s i m p l i c i t y of formulation since the l a t t e r required f i e l d data and estuary modeling to quantify i t s implications. 3.3 Nonlinear Programming Model This section describes the methodology of so l u t i o n of the NLP problem of section 3.1 by the combination of the simulation models of section 3.2 and a nonlinear optimization algorithm. 3.3.1 Methodology In section 3.2 the water q u a l i t y model was s p l i t into BOD and DO components r e s p e c t i v e l y . The solution methodology exploited t h i s separation and was as follows: 1) i d e n t i f y model parameters, boundary and i n i t i a l conditions; 2) simulate the BOD concentrations r e s u l t i n g from unit loadings i n a l l reaches over a l l time steps i n the decision horizon 3) use the BOD concentration (time and space) p r o f i l e s from 2) as input to a non-linear programming algorithm, with an embedded DO model, to i d e n t i f y the optimal BOD loadings subject to DO standards; 4) do 1) through 3) for s e n s i t i v i t y a n a l y s i s . Thus, steady (time invariant) unit BOD loadings were used i n step 2) above to simulate dynamic BOD p r o f i l e s which were l i n e a r l y combined i n step 3) as required for input into the NLP optimization algorithm. The objective function values were then determined for solutions i n the f e a s i b l e s p a c e d e f i n e d by the l e v e l o f c o m p l i a n c e e s t i m a t e d f r o m s i m u l a t e d DO p r o f i l e s , f o r a l l r e a c h e s o v e r a d e c i s i o n h o r i z o n . A f l o w c h a r t o f t h i s m e t h o d o l o g y i s p r e s e n t e d i n f i g u r e 3 .8 w h i c h combines a s p e c t s f o u n d i n t h e w o r k s o f G o r e l i c k (1982) and B i s h o p e t a l (1976) a p p l i e d t o a n o n l i n e a r programming p r o b l e m . 3 . 3 . 2 N o n - L i n e a r Programming A l g o r i t h m s T h e NLP p r o b l e m f o r m u l a t e d i n s e c t i o n 3 . 1 i s a c o n s t r a i n e d o p t i m i z a t i o n p r o b l e m w i t h l i n e a r o b j e c t i v e f u n c t i o n and n o n l i n e a r e q u a l i t y and i n e q u a l i t y c o n s t r a i n t s . A n e f f i c i e n t n o n l i n e a r o p t i m i z a t i o n a l g o r i t h m was t h e n s o u g h t . Due t o t h e c o m p u t a t i o n a l i n t e n s i t y o f the s o l u t i o n m e t h o d o l o g y p r o p o s e d i n s e c t i o n 3 . 3 . 1 , o n l y d e r i v a t i v e - f r e e m u l t i d i m e n s i o n a l o p t i m i z a t i o n a l g o r i t h m s were r e v i e w e d . A d d i t i o n a l l y , s i n c e a l l u n c o n s t r a i n e d methods c a n accommodate c o n t r a i n t s ( W a l s h , 1979) , the r e v i e w was r e s t r i c t e d t o s u c h m e t h o d s . H i m m e l b l a u (1972) s u r v e y e d some 15 d i f f e r e n t o p t i m i z a t i o n t e c h n i q u e s i n an o b j e c t i v e f ramework where a l l t e c h n i q u e s were e x e c u t e d on f i f t e e n t e s t p r o b l e m s and j u d g e d i n t e r m s o f : a) r o b u s t n e s s : w h e t h e r t h e a l g o r i t h m c o u l d s o l v e most o f t h e p r o b l e m s p o s e d t o w i t h i n a c e r t a i n p r e c i s i o n ; b) e f f e c t i v e n e s s : t h e number o f e v a l u a t i o n s o f t h e o b j e c t i v e f u n c t i o n i n a r r i v i n g a t t h e o p t i m a l s o l u t i o n ; c ) c o m p u t a t i o n t i m e : t h e t i m e t o t e r m i n a t i o n w i t h i n a d e s i r e d d e g r e e o f p r e c i s i o n . T h e d e r i v a t i v e - f r e e a l g o r i t h m s i n c l u d e d i n t h e t e s t ( w i t h t h e y e a r o f p u b l i c a t i o n ) w e r e : - 63 -INITIALIZATION (Rg 3.8o) SIMULATE BOD CONCENTRATIONS IN ALL REACHES DUE TO UNIT LOADINGS OVER THE DECISION HORIZON (Section 3 .13 .3 ; Fig 3.8b) I INITIAL GUESS OF BOD LOADINGS NONUNEAR PROGRAMMING ALGORITHM (Section 3.3.2) 1 LOADINGS OPTIMAL ? (CONVERGENCE CRITERIA MET ?) I OPTIMIZATION TECHNIQUE (Section 3.3.2) OBJECTIVE FUNCTION EVALUATION COMPUTE TOTAL BOD CONCENTRATION (Section 3.1.2; Fig 3.8c) SIMULATE DO CONCENTRATIONS (Section 3.2.3.4; Rg 3.I COMPUTE PROPORTION OF COMPLIANCE (Section 3.1.4; Fig 3.8e)  IF PROPORTION < ACCEPTABLE COMPUANCE SET PENALTY VALUE AS OBJECTIVE IF PROPORTION > ACCEPTABLE COMPUANCE COMPUTE OBJECTIVE VALUE (Eqn 3.1) Figure 3.8: Flowchart of the Solution Methodology for the NLP Programming Problem - 64 -FIGURE 3 . 8 a ) I N I T I A L I Z A T I O N I N I T I A L CONDITIONS : DISCHARGE, S T A G E , BOD, DO FOR A L L REACHES BOUNDARY CONDITIONS: UPSTREAM DISCHARGES, DOWNSTREAM T I D E S , UNIT BOD ADDITIONS PARAMETERS : BOD AND DO C O E F F I C I E N T S , LONGITUDINAL DISPERSION DATA : RIVER C R O S S - S E C T I O N S , BENTHIC DEMAND, PHOTOSYNTHESIS, TEMPERATURE FIGURE 3 . 8 b ) SIMULATE BOD CONCENTRATIONS SIMULATE yijt DUE TO A 200 kg/d UNIT LOADING OF BOD AS fijtCui.mg) I N FOUR REACHES SEPARATELY FOR 7 D A Y S , WHERE i IS THE ' I N P U T ' REACH AND j IS THE ' O U T P U T ' REACH ( i , j = 1 , . , 4 ) AND t IS FOR 10 MINUTE TIME PERIODS (SECTION 3 . 2 . 3 . 3 ) . REACH V E L O C I T I E S ARE DERIVED FROM THE HYbRODYNAMIC RIVER MODEL (SECTION 3 . 2 . 1 ) - 65 -FIGURE 3 . 8 c ) COMPUTE TOTAL BOD CONCENTRATIONS SIMULATE p j t AS 2 1 v i j t • x i / u i WHERE X i ARE THE BOD LOADINGS TESTED FOR OPTIMALITY I N THE NLP ALGORITHM AND u± IS 200 k g / d UNIT LOADING. FIGURE 3 . 8 d ) SIMULATE DO CONCENTRATIONS SIMULATE z j t AS gjtCPjt* 1^) WHERE gj t(.,.) I S THE DO MODEL DEVELOPED I N SECTION 3 . 2 . 3 . 4 . FIGURE 3 . 8 e ) PROPORTION OF COMPLIANCE ESTIMATE PROPORTION OF THE RATIO OF TOTAL Z j t ( j = 1 , . , 4 a n d t = 1 , . , 7 e v e r y 10 m i n ) WHICH EXCEEDS THE DO REGULATORY L I M I T (SECTION 3 . 1 . 2 ) . - 66 -1) Hooke-Jeeves (1961) 2) Nelder-Mead (1965) 3) Powell (1964) 4) Rosenbrock (1960) 5) Stewart (1967) The l a s t one was not pursued since i t required the computation of numerical d e r i v a t i v e s which would increase the computational burden. Based on computation time, the t e s t i n g procedure ranked the algorithms as: 1) Powell (6.1) 2) Hooke-Jeeves (9.6) 3) Nelder-Mead (10.0) 4) Rosenbruck (11.9) with t h e i r averaged ranking over a l l t e s t s r e l a t i v e to the other f i f t e e n algorithms tested, shown i n parentheses. The Powell algorithm appeared very promising for use i n t h i s study; however, due to the a v a i l a b i l i t y of computer code the Hooke-Jeeves algorithm was selected as the second best and implemented using FORTRAN code av a i l a b l e from Kuester and Mize (1973). 3.3.3 Hooke-Jeeves Method This d i r e c t search algorithm seeks the minimum of a mult i v a r i a b l e unconstrained nonlinear function F(x^1^2'''*xn^ which can be e a s i l y modified to incorporate constraints. The procedure assumes a unimodel function and thus should be used with multiple s t a r t i n g values i f t h i s assumption does not hold. The algorithm i s as follows (Kuester and Mize, 1973): - 67 -1) A base point i s picked and the objective function evaluated; 2) Local searches are made i n each d i r e c t i o n by stepping a distance s^ to each side and evaluating the objective function to see i f a lower function value i s obtained; 3) I f there i s no function decrease, the step size i s reduced and searches made from the previous best point; 4) I f the value of the objective function has decreased, a "temporary head" x. (k+1), i s located using the two i,o previous base points x^(k+1) and x^(k): x. (k+1) = x, (k+1) + a(x.(k+1) - x.(k)) l ,o J- i i where i i s the var i a b l e index = 1, 2, 3, N o denotes the temporary head k i s the stage index (a stage i s the end of N searches) a i s an acceleration f a c t o r , a >_ 1; 5) If the temporary head r e s u l t s i n a lower function value, a new l o c a l search i s performed about the temporary head, a new head i s located and the value of F checked. This expansion continues as long as F decreases; 6) If the temporary head does not r e s u l t i n a lower function value a search i s made from the previous best point; 7) The procedure terminates when the convergence c r i t e r i o n i s s a t i s f i e d . A flowchart of the above i s presented i n figure 3.9. Constraints are incorporated into the algorithm by i n d i c a t i n g a constraint v i o l a t i o n with a high F value (penalty) which would discourage exploratory and - 68 -Pick Base Point and Step Sizes. Evaluate Objective Function. Make Local Search Moving a Distance S. to Each Side and Evaluate Function Decrease Step Size No No Yes Locate Temporary Head ) > and Evaluate Function Yes Stop) Figure 3.7: Hooke and Jeeves (HOOKE ALGORITHM) Logic Diagram (Kuester and Mize, 1973) - 69 -p a t t e r n moves a c r o s s a c o n s t r a i n t b o u n d a r y ( W a l s h , 1 9 7 9 ) . The c o n s t r a i n t c a n be o f any c o m p l e x f o r m a s l o n g as i t c a n be a c c u r a t e l y d e t e r m i n e d . F o r t h e o p t i m a l BOD a l l o c a t i o n p r o b l e m t h e e v a l u a t i o n o f t h e j o i n t c o n s t r a i n t s was done as o u t l i n e d i n f i g u r e 3 . 7 a n d 3 . 8 e ) . The DO model was i n c o r p o r a t e d i n t o t h e H o o k e - J e e v e s a l g o r i t h m s i m u l a t i n g DO p r o f i l e s f o r a l l r e a c h e s a t e a c h o b j e c t i v e f u n c t i o n e v a l u a t i o n r e q u i r e d i n t h e p r o c e d u r e o u t l i n e d on t h e p r e v i o u s p a g e . Any n e g a t i v e Xj[ v a l u e s g e n e r a t e d i n t h e p r o c e d u r e were a l s o d i s c o u r a g e d w i t h p e n a l t y v a l u e s . T h e u n i m o d a l a s s u m p t i o n made by t h e H o o k e - J e e v e s method was t e s t e d w i t h the r e s u l t s p r e s e n t e d i n s e c t i o n 6 . 2 . T h e H o o k e - J e e v e s code was d e v e l o p e d and v a l i d a t e d w i t h sample p r o b l e m s u s e d by K u e s t e r and M i z e ( 1 9 7 3 ) . 3 . 4 Computer R e s o u r c e s and Code T h e m i c r o c o m p u t e r i s f a s t r e p l a c i n g m a i n - f r a m e m a c h i n e s i n t h e m o d e l i n g and o p t i m i z a t i o n o f s m a l l - t o m e d i u m - s i z e d e n g i n n e r i n g and management p r o b l e m s . A COMPAQ D e s k p r o 286 m i c r o c o m p u t e r (IBM P C - A T c o m p a t i b l e a t 8MHz) was u s e d t h r o u g h o u t t h i s s t u d y f o r the d e v e l o p m e n t and e x e c u t i o n o f computer code w r i t t e n i n FORTRAN 77 w i t h R y a n - M c F a r l a n d ' s RM/FORTRAN. Programs were l i m i t e d t o a 640 KB s i z e w h i l e t h e s c o p e o f t h e r i v e r m o d e l s was r e d u c e d t o 4 r e a c h e s ( f r o m an i n i t i a l 11) and a d e c i s i o n h o r i z o n o f one week ( r e d u c e d f r o m one month) s u c h t h a t e x e c u t i o n t i m e s were n o t p r o h i b i t i v e s i n c e t i m e s t e p - s i z e s o f 10 m i n u t e s were r e q u i r e d t o m a i n t a i n s t a b i l i t y o f t h e f i n i t e d i f f e r e n c e s o l u t i o n i n t h e NLP f o r m u l a t i o n o f s e c t i o n 3 . 1 . T h u s , a t o t a l number o f 4 r e a c h e s and 1008 t i m e p e r i o d s were u s e d . T h r o u g h o u t t h e work p e r f o r m e d , t h e m i c r o c o m p u t e r was e v a l u a t e d as t o i t s a p p l i c a b i l i t y t o t a c t i c a l a n d s t r a t e g i c management u s a g e i n t e r m s o f e x e c u t i o n t i m e and p r o b l e m s i z e . - 70 -4. FIELD INFORMATION 4.1 S i t e Description The a p p l i c a t i o n s i t e for t h i s study was the Nicomekl River i n the Nicomekl-Serpentine r i v e r basin i n the Municipality of Surrey i n B r i t i s h Columbia (see figu r e 4.1). The climate and geology was summarized from a 1963 'Drainage report of the Municipality of Surrey' presented below. 4.1.1 Climate "The climate of the southwestern section of the Municipality has been described as of the 'cool Mediteranean type, characterized by warm temperate rainy winters and cool dry summers'....while the remainder of the Municipality has 'a marine west-coast climate featured by warm temperate rainy winters and cool summers'. Surrey i s protected from the d i r e c t onslaught of storms from every d i r e c t i o n by various mountain ranges....the Coast Range to the north and the Cascades to the east. It i s protected from....rainstorms....from the west and southwest over the P a c i f i c Ocean by the mountains of Vancouver Island and the Olympic Mountains of northwest Washington. The e f f e c t of storms....from the southwest and west varies considerably....due to the u p l i f t caused by the Coast Mountains.... consequently the t o t a l annual r a i n f a l l and snowfall i s much greater i n the north-eastern section of the Municipality than i t i s i n the south-western section. The minimum to maximum average annual t o t a l p r e c i p i t a t i o n i s approximately 40 inches (1020 mm) to 70 inches (1780 mm)". The average annual temperature i s f a i r l y consistent throughout Surrey at approximately 50 degrees F (11.3 C e l s i u s ) . The ! F i g u r e 4-.lt L o c a t i o n p l a n s h o w i n g S u r r e y i n r e l a t i o n t o t h e G r e a t e r V a n c o u v e r a r e a ( D i s t r i c t o f S u r r e y , 1963) - 72 -average hours of sunshine vary from 1900 hours i n the southwest to 1500 hours i n the northwest section. "During the summer months of most years a l l the Municipality requires supplemental i r r i g a t i o n , p a r t i c u l a r l y during July to August, for the maximum development and y i e l d of crops." The Department of Agriculture and Food has suggested an i r r i g a t i o n duty of 18 inches (457 mm). "The worst drainage conditions occur when there are high r a i n f a l l i n t e n s i t i e s a f t e r a heavy f a l l of wet snow has f a l l e n on frozen ground which was previously saturated by r a i n " . Such an event was recorded i n 1935. "Wind e f f e c t s are not s i g n i f i c a n t i n terms of increasing evaporation from water surfaces but may caused extended period of high t i d e causing the r i v e r to back up". 4.1.2 Geology "Surrey i s part of the general area known as the Fraser Lowland which l i e s between two major geological features c a l l e d the Coastal Mountains and the Coastal Trough. The Surrey portion of the Fraser Lowland consists of extensive low h i l l s or uplands ranging up to 400 feet (122 m) i n elevation and separated by the f a i r l y wide f l a t bottomed v a l l e y s of the Nicomekl and Serpentine Rivers and the narrower v a l l e y of the Campbell River. The v a l l e y s of the former two r i v e r s were o r i g i n a l l y enlargements of the sea and were not caused by r i v e r erosion. The upper part of the Campbell River V a l l e y was probably created by the r i v e r as the land rose a f t e r being depressed r e l a t i v e to the sea by the l a s t i c e sheet.... between 5000 to 10000 years ago. It was estimated that the thickness of the i c e was 7500 feet (2286 m) causing great pressures - 73 -which account for the very hard nature of the clays and s i l t s which were l a i d down p r i o r to the i c e forming. "The upland areas of Surrey are of two main types. The f i r s t c onsists of a core of various types of clays, s i l t s , sands and gravels deposited by g l a c i a l i c e , streams, lakes, swamps and the sea. Thess uplands are usually f l a t or terraced with sufaces of g l a c i a l outwash found i n the Campbell River upland. The second type of upland has a s i m i l a r core but the areas are usually r o l l i n g and hummocky with surfaces of g l a c i a l t i l l and glacio-marine deposits occurring i n the Newton, Clayton and Sunnyside Uplands. "Bedrock l i e s several hundred feet below the above deposits and no bedrock i s to be seen anywhere i n Surrey. In the Boundary Bay area the bedrock i s over 2000 fee t (610 m) deep. The t i l l s , clays and s i l t s have low permeabilities discouraging high i n f i l t r a t i o n rates, The sand and gravel outwash allow higher i n f i l t r a t i o n ans sub-surface drainage. "The Nicomekl and Serpentine River Valleys consist mainly of peat underlain by impervious layers of clay and s i l t deposits. The groundwater body i s located near the surface and excess r a i n f a l l cannot be absorbed." - 74 -4.1.3 River Discharges, Temperature and I r r i g a t i o n The Nicomekl River i s a low l y i n g meandering t i d a l r i v e r as shown i n the photographs of fi g u r e 4.2. The r i v e r p r o f i l e i n fi g u r e 4.3 i d e n t i f i e s the break i n r i v e r slope i n the upland reaches. A d e t a i l e d plan of the watershed i n shown i n fi g u r e 4.4a) with measured cross-sections indicated. Summary data for cross-sections used i n t h i s study are shown i n f i g u r e 4.5b). H i s t o r i c a l monthly discharge data for the Nicomekl are available from Water Survey of Canada (Stn: 08MH105) and are summarized on an annual basis i n figu r e 4.5 showing 5%, mean and 95% discharge envelopes. A frequency analysis indicated that these monthly flows are log-normally d i s t r i b u t e d (figure 4.6). Stream temperatures from the Serpentine River for various stations and various sampling periods are summarized i n figure 4.7 by Ju l i a n day. These temperatures were applied to the Nicomekl i n l i e u of actual data to assess the temperature dependency of BOD and DO decay parameters. Following techniques used by several researchers surveyed i n sections 2.2.1 and 2.2.2, the r i v e r discharges and stream temperature could be generated by Monte Carlo or Markov techniques to account f o r temporal and possibly s p a c i a l v a r i a b i l i t y . This was not done i n t h i s study but would make an i n t e r e s t i n g extension. Detailed though sporadic water q u a l i t y data was av a i l a b l e for the Serpentine River i n a report by Moore (1985a). A s i m i l a r report for the Nicomekl i s i n preparation (Moore, 1985b). This grab sample data was from a monitoring network of stations on the r i v e r , although i t i s not s u f f i c i e n t f o r mathematical model development or v a l i d a t i o n . The use of t h i s data for model ' i d e n t i f i c a t i o n ' as discussed i n section 2.2.3 where - 75 -Figure 4 . 3 a : The Nicomekl River: Downstream Tidal Gate - 75a -Figure 4.3b: The Nicomekl River: Nicomekl at 40 Ave. - 75b -Figure 4 . 3 c : The Nicomekl River: Nicomekl before entering Mud Bay F i g u r e . %.4' N i c o m e k l R i v e r P r o f i l e N A • K Mitwr^. •«* «^  c A C*ft»5 . .- • N M R M p f « J e d f w * .• : \ , * '. ,-"f 'V. ' * •' *, 0MMLM MOM . B B t l ••••4. 0.6 i * " : ; * -• '-•" *eiu or itu* RieoMCKL s e n K N t m e " ftOOD CONTWH. M O J E t t ' tOC*TKW f l H F i g u r e 4.5*)*. P l a n of the N i c o m e k l - S e r p e n t i n e f l o o d p l a i n - 78 -4 CM O •—< X 0 -16 -12 -8 -A 0 A 8 12 16 20 E l e v a t i o n w i t h r e s p e c t to GSC ( f t ) F i g u r e 4^>feii)2 Graph of p e r i m e t e r v e r s u s r i v e r water e l e v a t i o n - 79 -N I C O M E K L R I V E R B E L O W M U R R A Y C R E E K MEAN MONTHLY D ISCHARGE E N V E L O P E S 300 2 B O + 2 0 O + S3 1 B O + o 1 0 0 + 8 0 + -+- -+- —i 1 — J A N F E B MAR A P R MAY J U N J U L A U G S E P O C T N O V D E C MONTHS LEGEND 5 % M E A N 9 5 % Figure 4.6 FREQUENCY OF LOG-NORMALIZED NICOMEKL M E A N M O N T H L Y F L O W S 0.30 0.06 H 0.04 • + 0.28 H 0.26 0.24 0.22 -0.20 -0.18 -0.16 0.14 -0.12 - + ffi 0.10 - D 0.08 -+ • • + + • °-02- t 5 • 0.00 -f ffl $ lP r ' 1 ' 1 ' 1 » 1 1 T * i -1 .350 -1.050 -0.750 -0.450 -0.150 0.150 0.450 0.750 1.050 Figure 4.7 L O G ( Q ? / Q a v g ) • D A T A + N O R M A L C U R V E SERPENTINE RIVER TEMPERATURES T 1 1 1 1 r 120 140 160 180 Figure 4.7: 1 1 1 1 1 1 1 1 1 1 1 1 D I 200 220 240 260 280 300 320 J U L I A N D A Y - 83 -model parameters and i n i t i a l conditions could be estimated, i s another i n t e r e s t i n g extension. DO data for the Nicomekl River was a v a i l a b l e for a l i m i t e d sampling i n 1974 to 1875 from Bourque and Hebert (1982) and was used to determine i n i t i a l and boundary conditions for the BOD/DO model. I r r i g a t i o n l i c e n s e information was a v a i l a b l e from the Water Management Branch giving maximum allowable seasonal ( A p r i l to September) i r r i g a t i o n uptakes. This provided an upper l i m i t of i r r i g a t i o n demands but was too coarse for t h i s study, even when monthly estimates were made (Van der Gulik, 1986). A d d i t i o n a l l y , actual uptake data was not available to characterize the frequency and volume d i s t r i b u t i o n s of withdrawals. Thus, both i r r i g a t i o n (and di t c h inflows) were not accounted f o r . This exclusion would be least v a l i d at low flow periods when i r r i g a t i o n withdrawals may have s i g n i f i c a n t impact on water q u a l i t y by reducing the d i l u t i o n and transport c a p a b i l i t y of the r i v e r . 4.1.4 Tides The tides i n Boundary and Mud Bays are mixed and mainly semi-diurnal d i s p l a y i n g a major M2 component with amplitudes of 0.8m and 156 to 164 phase (Crean, 1983). Sample t i d e records are shown i n figure 4.8. 4.1.5 Land Use and P r i n c i p l e J u r i s d i c t i o n s "The Serpentine-Nicomekl lowlands have been reserved f o r a g r i c u l t u r e and the (Surrey) Council has adopted a p o l i c y aimed at protecting these farmlands against frequent flooding. On the Upland areas, the r u r a l and forested lands are s t e a d i l y being occupied by housing and commercial development." (Figure 4.9) TIDAL GAUGE HEIGHT AT TSAWASSEN A U G U S T 1 9 7 7 ( f t - G S C D A T U M ) 0 4 8 9 6 1 4 4 1 9 2 2 4 0 2 8 8 3 3 6 3 8 4 4 3 2 4 8 0 5 2 8 5 7 6 6 2 4 6 7 2 7 2 0 Figure 4 . 9 H O U R S - 85 -g u r e 4.10 P l a n s h o w i n g e x p e c t e d u l t i m a t e l a n d u ( D i s t r i c t o f S u r r e y , 1 9 6 3 ) - 86 -The a g r i c u l t u r a l lands are are protected by a dyking system maintained by the New Westminister Dyking Commission. At the downstream ends of both r i v e r s are control structures i n the form of t i d a l gates which minimize s a l i n i t y i n t r u s i o n and attenuate t i d a l f l u c t u a t i o n s . The present structures were i n s t a l l e d by the federal Ministry of the Environment i n 1972. A seasonal i r r i g a t i o n l i c e n s i n g system for the r i v e r waters i s maintained by the Water Management Branch of the p r o v i n c i a l Environment Ministry. This agency grants licenses and has the mandate to control actual i r r i g a t i o n (by l i c e n s e s e n i o r i t y ) which may occur towards the end of summer when i r r i g a t i o n demands are high and r i v e r flows low. The t i d a l f l a t s area of the watershed i s an e c o l o g i c a l sanctuary and i s supported by the p r o v i n c i a l Fish and W i l d l i f e Branch. Both the Nicomekl and the Serpentine r i v e r s support salmon runs which t r a v e l through Boundary Bay to spawn. The federal F i s h e r i e s and Oceans i s resposible for protecting the f i s h e r y resource and i n conjunction with the p r o v i n c i a l Water Quality Branch i s pursuing strategies to a l l e v i a t e and control the negative impact of poor water q u a l i t y . It i s the mandate of the Municipality of Surrey to formulate land use p o l i c i e s while providing basic services and i n f r a s t r u c t u r e to the residents, businesses and i n d u s t r i e s . The a g r i c u l t u r a l land i s also part of the p r o v i n c i a l A g r i c u l t u r a l Land Reserve with land sales and purchase regulated by the Land Use Commission. The p r o v i n c i a l Ministry of Agriculture and Food also a s s i s t s a g r i c u l t u r a l land use decisions by providing technical support and consultation (Van der Gulik, 1983). The primary source of poor water q u a l i t y i n the Serpentine-Nicomekl r i v e r s i s the long and short term non-point organic waste p o l l u t i o n from a g r i c u l t u r a l lands bordering the r i v e r s and secondary - 87 -short-term commercial and i n d u s t r i a l point sources. Currently the Serpentine r i v e r has the poorer water q u a l i t y where s i g n i f i c a n t f i s h k i l l s have been reported as a r e s u l t of low dissolved oxygen (Moore, 1985a). The Nicomekl r i v e r i s estimated to reach t h i s same condition i n 'about f i v e years' (Gough, 1985). These organic wastes appear to reach the the r i v e r s p r i m a r i l y by surface runoff processes and groundwater drainage with groundwater seepage playing a secondary r o l e . The low dissolved oxygen problem has reach c r i s i s proportion i n the Serpentine r i v e r such that a f e d e r a l l y funded project i s underway to aerate a c r i t i c a l lowland r i v e r reach. (Mavinic, 1986). The Nicomekl which was chosen as the r i v e r of i n t e r e s t f o r t h i s study, due to r i v e r data a v a i l a b i l i t y , though a p p l i c a t i o n of the methodology presented previously may have had more immediate p o l i c y implications i f applied to the Serpentine. 4.5 Data Summary The following data sets were av a i l a b l e for t h i s study: a) monthly and annual Nicomekl r i v e r flows for the purpose of quantifying the longer-term hydrologic regime; b) hourly t i d a l data from the Tsawassen t i d a l s t a t i o n number 2/4900/12307 f o r 1972 and 1977; c) hourly stream flow data from Nicomekl r i v e r stations WSC-08MH105 and WSC-08MH050 below Murray Creek and at 192nd Street r e s p e c t i v e l y for 1972 and 1977; d) cross-sectional data from various sections on the Nicomekl River from 1970 Environment Canada project data; - 88 -e) hourly t i d a l predictions for Tsawassen for twenty years provided by the I n s t i t u t e of Ocean Sciences at P a t r i c i a Bay, B r i t i s h Columbia. f) stream water q u a l i t y measurements for the Serpentine r i v e r from monitoring surveys done by the Water Quality Branch and the Environmental Protection Service. g) i r r i g a t i o n l i c e n s e information on the maximum seasonal withdrawals allowed from the Water Management Branch. h) water q u a l i t y information for the purposes of v a l i d a t i n g the water q u a l i t y model developed i n section 3.2.2 was a v a i l a b l e i n the , report of Pence et a l (1968) for the Delaware Estuary. Data sets b), c) and d) were used to c a l i b r a t e and v a l i d a t e the water quantity model of section 3.2.1. Data set f) was used to i d e n t i f y c r i t i c a l low dissolved oxygen periods and for stream temperature data. Since s u f f i c i e n t water q u a l i t y data was not av a i l a b l e to v a l i d a t e the BOD/DO model data set h) was used as one method of p a r t i a l v a l i d a t i o n while model parameters were obtained from the l i t e r a t u r e . - 89 -5. MODEL CALIBRATION AND VALIDATION The models developed i n section 3.2 were c a l i b r a t e d and vali d a t e d using data sets described i n section 4.2. This chapter summarizes the v a l i d a t i o n procedure j u s t i f y i n g the use of the models, although v a l i d a t i o n of the water q u a l i t y model on another r i v e r does not va l i d a t e the model for use on the s i t e of i n t e r e s t (Ambrose, 1986). Transient conditions were used for v a l i d a t i o n so as to ensure a p p l i c a b i l i t y to the dynamic regime i n question. 5.1 Hydrodynamic Model C a l i b r a t i o n runs were made using the data sets described i n section 4.2a) and b). S p e c i f i c a l l y , the boundary condtions were established by: a) t i d a l gauge height data from Tsawassen; b) discharge data for the st a t i o n below Murray Creek which was computed from a nonlinear regression of a stage-discharge curve (Jamal, 1980); c) i n i t i a l conditions were taken as steady state flow of 50 c f s and 0.2 f t depth of flow as a s i m p l i f i e d assumption. The transient e f f e c t of the i n i t i a l conditions were allowed to di s s i p a t e over the f i r s t 2000 i t e r a t i o n s of the model. Nicomekl cross-section data of section 4.2d) was used to characterize the r i v e r geometry. The Manning's 'n' was used as the c a l i b a r a t i n g parameter and manipulated during c a l i b r a t i o n u n t i l a 'good' v i s u a l f i t was obtained with A p r i l 1977 data. This was compared to f i e l d estimated values of s i m i l a r r i v e r s of 0.03 to 0.05 (Chow, 1959). - 90 -The d i s c r e p a n c i e s between t h e s i m u l a t e d and a c t u a l r i v e r s t a g e s may be due t o : a ) not a c c o u n t i n g f o r s u r f a c e w a t e r d i s c h a r g e f r o m d i t c h e s and m i n o r t r i b u t a r i e s and i r r i g a t i o n u p t a k e ; b) g r o u n d w a t e r r e c h a r g e and d i s c h a r g e a t d i f f e r e n t p o i n t s i n t h e t i d a l c y c l e d e p e n d i n g on t i d a l v e l o c i t y and s o i l m o i s t u r e c o n d i t i o n s a t d i f f e r e n t l o c a t i o n s a l o n g t h e r i v e r . M o d e l v a l i d a t i o n was made f o r A p r i l 1972 as shown i n f i g u r e 5 .2 and c o n f i r m e d by s t a t i s t i c a l t e s t i n g ( t a n d F t e s t s a t t h e 95% c o n f i d e n c e l e v e l ) . An i s s u e n o t a d d r e s s e d h e r e i s t h e v a r i a b i l i t y o f t h e c a l i b r a t i n g p a r a m e t e r w h i c h c o u l d be d e t e r m i n e d , b o t h f o r i n t r i n s i c a n d s e a s o n a l v a r i a t i o n , u s i n g r e a l - t i m e e s t i m a t i o n t e c h n i q u e s as o u t l i n e d i n s e c t i o n 2 . 2 . 4 . C a l i b r a t i o n o f t h e m o d e l was done f o r A p r i l 1977 d a t a f o r h i g h e r d i s c h a r g e s s t h a n summer f l o w s b u t where t h e d y n a m i c s o f f l o o d h y d r o g r a p h s w o u l d be a good t e s t o f t h e m o d e l . Summer f l o w s may i n d e e d be b e t t e r f i t w i t h d i f f e r e n t v a l u e s o f ' n ' . I n summary, t h e h y d r o d y n a m i c m o d e l a p p e a r s t o be a ' w e l l i d e n t i f i e d ' m o d e l . 5.2 BOD/DO M o d e l T h e v a l i d a t i o n o f t h e BOD/DO m o d e l f o r t h e N i c o m e k l R i v e r was n o t p o s s i b l e due t o t h e l a c k o f s u i t a b l e f i e l d d a t a . A l t e r n a t e methods f o r v a l i d a t i o n o f t h e m o d e l were r e v i e w e d : a) v a l i d a t i o n o f t h e a l g o r i t h m by s i m i l a r s t u d i e s documented i n the l i t e r a t u r e ; b ) v a l i d a t i o n o f t h e a l g o r i t h m b y c o m p a r i s o n w i t h w e l l - k n o w n a n a l y t i c m o d e l s ; c ) v a l i d a t i o n o f t h e m o d e l on a d a t a s e t f r o m a n o t h e r r i v e r s y s t e m . VALIDATION OF THE HYDRODYNAMIC MODEL 2 3 7 1 1 1 9 1 6 7 2 1 5 2 6 3 3 1 1 3 5 9 4 0 7 4 5 5 5 0 3 5 5 1 5 9 9 6 4 7 6 9 5 Figure 5.2 HOURS SIMULATED X RECORDED - 92 -5.2.1 V a l i d a t i o n i n the L i t e r a t u r e A review of studies using the model described i n section 3.2.2 revealed a few studies which have va l i d a t e d the use of t h i s model for f i e l d conditions. B e l l a (1970) used the model (in two-dimensions) to study the DO p r o f i l e s i n Lake Sammamish, Washington. The processes of reaeration, photosynthetic oxygenation, v e r t i c a l mixing and oxygen uptake were modeled while accounting for temperature v a r i a t i o n . Figure 5.2 shows the r e s u l t s from t h i s study with 5.2b) and 5.2c) as the measured and predicted v e r t i c a l DO p r o f i l e at two lake s t a t i o n s . S e n s i t i v i t y analysis was also performed for the c o e f f i c i e n t of atmospheric reaeration, photosynthetic and r e s p i r a t i v e rates and the v e r t i c a l dispersion c o e f f i c i e n t . Aiba and Ohtake (1977) used the same modeling procedure i n t h e i r study of the balance of P04-P (phosphate) i n the shallow p o l l u t e d Tama-gawa River, Japan. They modeled the processes of convection, dipersion, suspended and river-bed adsorption, hydrolysis, benthic demand, aquatic demand and decomposition of r i v e r organics. This sophisticated material balance model was used to predist the P04-P concentrations at four r i v e r stations as validated i n f i g u r e 5.3. These two studies address d i f f e r e n t p o l l u t i o n problems s u c c e s s f u l l y when using the basic B e l l a and Dobbins (1968) model. Though these studies do not v a l i d a t e the BOD/DO model for the r i v e r of concern i n t h i s study, they do lend c r e d i t and j u s t i f i c a t i o n to the use of the algorithms proposed. - 93 -FIG. 53,—DO PROFILES IN RESPONSE TO DIFFERENT PARAMETERS (Bella, 1970) 0.1 0.2 sT 0.0 9 11:00 ° 12:50 11:30 16 :00© 10.0 35.0 30.0 25.0 20.0 DISTANCE FROM ESTUWY (WO Fig.54. P 0 4 - P concentrations, observed (•) and simulated ( ° ) . (August 18, 1972); time of sampling is shown in the figure. (Alba and Ohtake, 1977) - 94 -5.2.2 V a l i d a t i o n Against A n a l y t i c Models A s i m i l a r v a l i d a t i o n process as presented by B e l l a and Dobbins (1968) was developed where two scenarios were considered: i ) slug BOD r i v e r loading of a t i d a l r i v e r to t e s t the convection, dispersion and decay processes of the BOD model; i i ) continuous BOD estuary loadings to tes t steady-state BOD and DO p r o f i l e s . For both scenarios, r e s u l t s using the model developed i n section 3.2.2 were compared with well known an a l y t i c models discussed by B e l l a and Dobbins (1968) and t h e i r own numerical r e s u l t s . 5.2.2.1 Slug Loads An i n i t i a l load of 100 un i t s i n a sinusoidal t i d a l r i v e r was used with parameters summarized i n table I Table I : PARAMETERS USED FOR ANALYTICAL VALIDATION Parameter Value Unit Remarks •2—r_ -1 D 2 mi .day Dispersion k 1 0.3 day" BOD Decay A x n 0.25 miles Reach Length A T 1/96 days Time Step U 17.3sin(2tt;t/P) mi.day" T i d a l V e l o c i t y P 12.5 hours T i d a l Period a miles Slug Length An a n a l y t i c a l s o l u t i o n f o r the instantaneous l i n e source heat conduction problem was presented by Carslaw and Jaeger (1960, as reported by B e l l a and Dobbins, 1968; and Sageev, 1986) as: L(x,T) = 50 [ erf(a-(x-x')) - erf(a+(x-x')) ] e ~ k l T 2 DT 2 DT where x' = j U(t) dt. 'o A favourable comparison with B e l l a and Dobbins' (1968) r e s u l t s and the a n a l y t i c a l model are shown i n table II . - 95 -Table | 1 : SLUG LOAD COMPARISONS Distance 1 T i d a l Cycle 2 T i d a l Cycles 6 T i d a l Cycles (mi) A* B+D** s*** A B+D S A B+D S 0 5.903 5.900 5.902 3. 572 3. 570 3.572 1. 104 1. 103 1. 104 1 4.646 4.647 4.648 3. 169 3. 168 3. 169 1.061 1.060 1.061 2 2.266 2.268 2.268 2. 212 2. 212 2.213 0.941 0.940 0.941 3 0.684 0.683 0.683 1. 215 1. 215 1.216 0.770 0.771 0.771 5 0.015 0.014 0.014 0. 179 0. 177 0.178 0.406 0.406 0.407 7 0.000 0.000 0.000 0. 010 0. 010 0.010 0. 155 0. 156 0. 156 * A n a l y t i c a l Model ** B e l l a and Dobbins (1968) *** This Study. For steady-state BOD-DO p r o f i l e s , O'Connor developed the following a n a l y t i c a l equations: L(>0 - L© €*P [_* X ( 5 ) * ] > U = ^  J COO .Q - Fuf e r f ± X ( § ) * ] « e [ i X ( f )*] j , F . K J ' " 3^— 1% I These were used for v a l i d a t i o n of the steady DO and BOD concentrations due to steady loadings. Table Vl| shows a favourable comparison for BOD and DO for steady state a f t e r 19 t i d a l cycles. Longer simulations would have yielded c l o s e r r e s u l t s . Table I | | : STEADY STATE COMPARISONS stance BOD (mg/1) DO > (mg/1) (mi) A B+D S A B+D S 1 6. 653 6.671 6.674 1. 833 1.714 1. 724 2 4. 516 4.522 4.530 2. 538 2.420 2. 427 3 3. 066 3.065 3.074 3. 373 3.255 3. 262 4 2. 082 2.077 2.086 4. 200 4.082 4. 088 6 0. 959 0.952 0.961 5. 592 5.479 5. 487 10 0. 204 0.200 0.204 7. 175 7.087 7. 096 16 0. 020 0.018 0.020 7. 872 7.831 7. 837 The two comparisons above resulted i n maximum deviations from the analytical solution of 0.6% and 5.9% respectively. - 96 -The comparisons above v a l i d a t e the code developed i n t h i s study as being equivalent to the a n a l y t i c slug load equations and that of B e l l a and Dobbins (1968). 5.2.3 V a l i d a t i o n Using Alternate F i e l d Data Data for the Delaware Estuary was a v a i l a b l e from the study of Pence et a l (1968), who used c h a r a c t e r i s t i c methods i n solving a f i n i t e - d i f f e r e n c e scheme of BOD and DO dynamics i n the estuary. This data was used i n an attempt to demonstrate the a p p l i c a b i l i t y and v a l i d i t y of the model developed i n section 3.2.3. A sucessful v a l i d a t i o n here does not imply that the BOD/DO i s validated for the Nicomekl River. A f i e l d data gathering program based on a r a t i o n a l i z e d experimental design for the Nicomekl would be required to throroughly v a l i d a t e the model, or more appropriately, to " i d e n t i f y " the best of various models av a i l a b l e using estimation techniques. Pence et a l (1968) present the following data i n t h e i r study: a) d a i l y r i v e r flow data at a Delaware gauging s t a t i o n ; b) d a i l y r i v e r temperature at another s t a t i o n ; c) mid-tidal r i v e r cross-section data for 30 sections; d) steady-state BOD loadings (carbonaceous and nitrogenous) i n l each reach; e) equations describing the temperature dependency of BOD decay and reaeration c o e f f i c i e n t s ; f) steady benthic oxygen demand for each reach (photosynthesis was n e g l i g i b l e ) ; and g) base (20 deg. c e l c i u s ) decay and dispersion c o e f f i c i e n t s . I n i t i a l conditions were not s p e c i f i e d while boundary conditions were - 97 -established by the r i v e r discharge. No diurnal t i d a l data was a v a i l a b l e so that mid-tidal cross-sectional data was used. There were no c a l i b r a t i n g parameters. Using data a) through g), DO concentrations i n the Delaware for 1964 were simulated with time steps of 0.25 days. The three v a l i d a t i o n cross-sections used by Pence et a l (1968) are shown i n figures 5.4 a ) , b), c ) , with recorded data, model r e s u l t s of these researchers and the simulation performed i n t h i s study. It i s apparent that the BOD/DO model of section 3.2.3 predicts more conservative p r o f i l e s . This i s l i k e l y due to the averaged t i d a l conditions used ( e s p e c i a l l y f or section 22 which the furthest downstream). The upstream section 4 and 7 y i e l d better f i t s to the data as the t i d a l influence i s reduced. The model r e f l e c t s a greater s e n s i t i v i t y to r i v e r temperature than the Pence et a l (1968) model. This v a l i d a t i o n indicates the usefulness of the BOD/DO models developed i n t h i s study, although f u l l v a l i d a t i o n was not obtained. As mentioned previously, v a l i d a t i o n of water q u a l i t y models with l i m i t e d i l l - d e s i g n e d data programs i s a d i f f i c u l t excercise. 5.2.4 Summary In summary, the BOD/DO model has not been f u l l y v a lidated f or the r i v e r system under study since: a) the a l t e r n a t i v e data set and ap p l i c a t i o n studies surveyed dealt with other r i v e r systems; b) the time scales of the phenomena analyzed by the studies i n a) are of d i f f e r e n t orders: seasonal versus d i u r n a l ; c) the a n a l y t i c models only v a l i d a t e the algorithm used i n the BOD/DO model but do not v a l i d a t e i t s a p p l i c a b i l i t y . DISSOLVED OXYGEN PROFILE AT DELAWARE ...MEMORIAL BRIDGE (SECTION 22) c ft w o 1 H H-M* 01 ' a. oi 00 o n t—* . «s re oo re 3 re CD c in 1 12 11 ^ 7 e W 6 O X O 5 Q u > g 4 / 1 7 [ 1 \ v.. / / / / # N / / » *fc> \ \ / • • • \ \ \ ^ — THIS STUDY • V • > \ \ \ \ • * \ \ OR • / i i < tt • X \ j— BOAT RUN • • [ / ' I r I 1 /. • r • \* • 1. \ 1 T • / i i i • \ • \^  DMPUTED i 1 . t VO 00 (Pence et a l , 1968) . JAN. FEB. MARCH, , APRU. , MAY j, JUNE : JULY: . AUG. „ SEP.. , , OCT. NOV; DEC. DISSOLVED OXYGEN PROFILE AT TORRESDALE (SECT. 7) VO (Pence et a l , 1968) JAN. FEB. MARCH APRIL MAY JUNE JULY AUG. SEP. OCT. NOV. DEC. 1964 F i g u r e 5.5b". D i s s o l v e d Oxygen a t T o r r e s d a l e , Pa. PROGRAM VALIDATION DISSOLVED OXYGEN PROFILE AT BRISTOL (SECTION 4) - 101 -d) there i s i n s u f f i c i e n t f i e l d water q u a l i t y data to v a l i d a t e the BOD/DO model for the Nicomekl r i v e r . The a p p l i c a b i l i t y and use of the model must therefore be used with due caution and the r e s u l t s presented i n the next chapter, taken as estimated predictions subject to f u l l model v a l i d a t i o n . - 102 -6. MODEL APPLICATION TO THE NICOMEKL RIVER This chapter discusses the a p p l i c a t i o n of the models developed i n chapter 3, following the methodology of figure 3.8, to the study area i n solving the problem formulated i n section 3.1. S p e c i f i c a l l y the methodology used was as follows: 1) choose model parameters for model a p p l i c a t i o n (section 6.1.5) with discharges at the 5% and 95% percentage l e v e l s ; 2) using the BOD model (section 3.2.2.3) and the hydrodynamic model (section 3.2.1) simulate, with models parameters and discharges above, BOD concentration time p r o f i l e s r e s u l t i n g from unit BOD loadings, for a l l reaches i n the r i v e r for each time period i n the decision horizon; 3) using the NLP algorithm (section 3.3.3) and DO model (section 3.2.2.4) obtain the optimal BOD loadings which y i e l d DO p r o f i l e s which meet the regulatory l i m i t s at some r i s k l e v e l ; 4) test the s e n s i t i v i t y of the optimal solutions to model parameters (section 6.3). 6.1 Factors i n BOD Simulation A v a r i e t y of factors were considered i n the p r e d i c t i o n of BOD p r o f i l e s . 6.1.1 River Discharge Steady and unsteady discharges were analysed. Using data f o r 1971 and 1977, the i n t e r a r r i v a l times between storm events were estimated and summarized i n a frequency p l o t of fi g u r e 6.1. Figure 6.1: Storm Interarrival Time for the Nicomekl River below Murray Creek for 1972 and 1977 data - 104 -The mean i n t e r a r r i v a l time of about 67 hours (2.8 days; i n a log-normal d i s t r i b u t i o n ) was about f o r t y percent of the one week decision horizon. Thus the use of a synthetic hydrograph generator developed i n t h i s study did not appear to be j u s t i f i e d for the purposes of simulating upstream inflows since the v a r i a b i l i t y of the synthetic hydrographs over one week would require extensive use of Monte Carlo simulation techniques extending considerably the execution times of 'production' runs of the model beyond the time l i m i t s for completion of t h i s study. A d d i t i o n a l l y , since c r i t i c a l periods occur during the dryer summer months, the use of steady discharges appeared s u f f i c i e n t to describe t h i s boundary condition. These discharges could be modeled either as constants or as some stationary autocorrelated process. The former s i m p l i f i e d method was used though the l a t t e r could have been e a s i l y incorporated as was done by Jamal (1980). 6.1.2 Tides A one week t i d a l period was chosen to represent t y p i c a l conditions at the downstream boundary of the r i v e r system. This was required since d i f f e r e n t t i d a l periods may r e f l e c t d i f f e r e n t t i d a l harmonics of long or short frequencies. The use of various t i d a l periods could have been made to r e f l e c t maximum and minimum t i d a l conditions as s e n s i t i v i t y tests on the optimal decisions. This was not done i n t h i s study, again, due to lengthy model execution times required on the microcomputer. An a l t e r n a t i v e source of t i d a l data f o r the downstream boundary was that a v a i l a b l e from the I n s t i t u t e of Ocean Sciences. This hourly data was generated by a model based on Gordon's t i d a l heights analysis by using harmonic constituents for s t a t i o n of i n t e r e s t . Predictions f o r - 105 -Tsawassen were made using 62 constituents including a high astronomical frequency of 18.6 years. This data was i n i t i a l l y intended for use i n a Monte Carlo simulation framework which was l a t e r discarded due to constraints on time and the computing resources a v a i l a b l e . Data was also a v a i l a b l e from coastal models developed by Crean (1983) which simulated the mixed tides of the S t r a i t s of Georgia and Juan de Fuca though the primary focus of the study was to model estuarine c i r c u l a t i o n . 6.1.3 BOD Loadings Figure -6.la shows the r i v e r reaches considered f o r BOD loading which were chosen to correspond to the scheme used for the hydrodynamic model. Continuous unit BOD loadings of 200 kg/day were assumed for each of the four reaches although stochastic loadings of the same mean could have been assumed. An i n i t i a l choice of eleven reaches for consideration was reduced to four due to the necessarily lengthy microcomputer execution times. A s t a b i l i z a t i o n period of 24 simulated hours was used for both the hydrodynamic and the BOD/DO models to d i s s i p a t e the e f f e c t of the i n i t i a l conditions. From the v a l i d a t i o n s i n chapter 5, t h i s period appeared to be adequate. Decay and l o n g i t u d i n a l dispersion c o e f f i c i e n t s could also be considered as stochastic, and generated using means and variances from appropriate s t a t i s t i c a l d i s t r i b u t i o n s and including possible autocorrelations. This was not done i n t h i s study. 6.1.4 Photosynthesis Data from the l i t e r a t u r e was used to determine an oxygen source l e v e l consistent with the observations of Moore (1985a) and Bourque and Hebert (1982). Steady source rates were assumed as s i m p l i f y i n g F i g u r e 6 . 1 a : . BOD/DO M o d e l Schemat ic o f R i v e r R e a c h e s B O P / D O M O D E L PREDICTIONS : REACH 1 REAC14 2 REACH 3 REACH 4 h< >^ >+« X DISTANCE. (-ft) 0 8 0 0 1600 7400 13200 35800 3 O 0 0 3 8 0 0 ACXaCO 43800 Al&yO 5)810 0 P 5 T P E A M | I I I I I I I I I I I o O N M T P P O P Y H A M I C M O D E L P p E D I C T l O H S •  U Q 14 <3 I I I L C«OS&-seOTIOH: 1 2. I i T l X A L TIPAL dCUHVm GATE Q J_ 14 Q Q 6 10 14 _L_ l i 1Z 1 W5C STATION Q = DISCHARGE H =• HT. OF PUCW BELOW MURRAY CKEEK - 107 -assumptions although more complex models for diurnal v a r i a t i o n s are ava i l a b l e (Jorgensen, 1983). An equivalent rate of 0.09 mg/l/day was used for each r i v e r reach 6.1.5 DO Regulatory Limit The regulatory l i m i t , c r e g » depends on the specie of f i s h which frequent the r i v e r for spawning and feeding. A review of minimum DO requirements for f i s h species by Davies (1975) i d e n t i f i e d a l e v e l of 5.0 mg/1 for the highest protection l e v e l 'A' for salmonids. This l e v e l of protection allows one standard deviation from the 'mean average i n c i p i e n t oxygen response threshold". A more exacting protection regime would also account for ranges i n r i v e r v e l o c i t y and minimum depths of flow required, say, during spawning. These standards would be more relevant to ad d i t i o n a l studies of flow augmentation and i r r i g a t i o n where multivariable controls were assessed. 6.1.6 Parameters Used The parameters used for the simulation and optimization runs made are summarized below: MODEL PARAMETER VALUE UNITS REMARKS Decision Horizon Hydrodynamic Water q u a l i t y n A t 7 day reduce from 30 days due to computational burdens. 0.060 - from c a l i b r a t i o n 30 sec from s t a b i l i t y c r i t e r i a 50 c f s i n i t i a l condition 0.2 f t i n i t i a l condition 2 1/day Dobbins(1964), Orlob(1983), Ott (1976) - 108 -k 2 1/day Orlob(1983), 1/2 3/2 [Du ' /H ' ] * * Ott(1976) D 2 sq mi/d Dobbins(1964) D b 0.5 gm/d/m2 Ambrose(1986) A T 600 sec s t a b i l i l t y c r i t e r i a BOD Q/BOD^* 0 mg/1 boundary conditions D O q / DCh* 11.5 mg/1 boundary conditions C 5.0 mg/1 Davies (1975) reg NLP model obj. function c o e f f i c i e n t s -for: X1 11600 f t length of reach 1 X2 17000 f t length of reach 2 X3 10000 f t length of reach 3 X4 7250 f t length of reach 4 * non-pristine background conditions for BOD and DO could have been e a s i l y incorporated. ** This i s a reaeration model proposed by O'Connor (1958) where D i s -5 2 the oxygen d i f f u s i v i t y through l i q u i d f i l m (8.1x10 f t /hr) and u i s the mean t i d a l v e l o c i t y and H i s the mean depth. 6.2 Test of the Unimodal Assumption The key assumptions made i n the Hooke-Jeeves algorithm i s the unimodality and smoothness of the constraints which bound the f e a s i b l e region. This ensures the l o c a t i o n of a glo b a l optimal. For multimodal constraint surfaces only l o c a l optima can be guaranteed, r e q u i r i n g the use of multiple s t a r t - p o i n t s (Kuester and Mize, 1973). In t h i s study, the constraints are nonlinear and need exploration to determine unimodality. For a constant discharge of 25.0 cfs and a regulation l i m i t of 5.0 mg/1 (k^ = 0.3 per day and D= 2 mi^ per day) the f e a s i b l e s o l u t i o n space for t h i s problem was explored by - 109 -g r i d search. The search was easier to v i s u a l i z e due to the manageable number of variables (four) which would not have been the case for a higher dimensioned problem. For various l e v e l s of loadings i n reaches 1 and 2 the estimated r i s k l e v e l was determined for two loading l e v e l s i n reach 3. The loading i n reach.4 appeared to consistently v i o l a t e the constraints except at n e g l i g i b l e load l e v e l s and was assumed to be zero. This l i m i t e d search indicated that the constraints bounding the so l u t i o n space were apparently unimodal. However, there were portions of the constraint surface which were almost p a r a l l e l to the objective function, posing possible convergence problems. This was confirmed with further t e s t i n g such that a reduction i n convergence c r i t e r i a from 1% to 0.001% yielded better optimal points. Figure 6.3 shows the constraints on the s o l u t i o n space for various compliance l e v e l s for BOD loadings i n reaches 1, 2 and 3 (from cross-section 3 to 5, 5 to 7 and 7 to 9 respectively) r e s u l t i n g from a f u l l e r g r i d search. The constraints indicated possible l o c a l multi-modality suggesting the use of multiple s t a r t points. The objective function i s also shown being almost p a r a l l e l to the constraints near the optimal region of the s o l u t i o n space. 6.3 Results and S e n s i t i v i t y Analysis The optimization problem formulated was solved using the methodology previously o u t l i n e d . The s e n s i t i v i t y of the optimal decision to v a r i a t i o n s i n the model parameters were tested as shown i n table IV • - 110 -0 ZOO 400 600 800 IOO0 I20O I400 1600 IflOO 2000 X 2 (Htfl Figure 6.3: Graphical Representation of a Limited Grid Search of the Feasible Region Showing Various Compliance Constraint Surfaces - 111 -Table "Iv*: S e n s i t i v i t y Ranges for Parameters Tested. PARAMETER VALUES UNITS da y - 1 0.05, 0.3 0.5 , 2.0 mi /day f t J / s e c mg/1 Q C iance Level 5, 25, 37.5 3, 5, 7, .8, .9, 1.0 Multiple s t a r t points were used to ensure consistency i n the optimal solutions and detection of small-scale anomalies i n the s o l u t i o n surface possible y i e l d i n g l o c a l optima. The NLP algorithm took an average of 1 hour when executed on the COMPAQ 286, to^generate an optimal s o l u t i o n . With f i v e sets of multiple s t a r t points and convergence c r i t e r i a , optimal BOD loadings were obtained. These may s t i l l , however, r e f l e c t l o c a l optima, re q u i r i n g more i t e r a t i o n s of the NLP algorithm, to help i d e n t i f y the global optimal s o l u t i o n . The NLP algorithm may also be unstable near the optimal point where the constraint boundary and the objective function are almost p a r a l l e l . This problem-specific c h a r a c t e r i s t i c may well be algorithm independent (of the ones reviewed) and may require incresed number of s t a r t points for adequate r e s o l u t i o n . However, tentative conclusions were drawn from computer runs for discharge l e v e l s of 5, 25 and 37.5 c f s , c r e g °f 3, 5, and 7 mg/1 and compliance l e v e l s of 0.8, 0.9 and 1.0. 6.3.1 Discharge Levels the three steady discharge l e v e l s used (5, 25 and 37.5 c f s ) . The trend i n BOD loadings appears consistent with i n t u i t i o n , namely higher loadings allowable at higher discharges due to added d i l u t i o n . In going Figures 6.3 to 6.5 present the objective function values for - 112 -Sensitivity of Optimal Loadings to Discharge Levels 0 1.0 COMPUANCE A 0-9 COMPLIANCE. V 0-6 COMPLIANCE REACH 1 REACH 2 REACH 3 OBJECTIVE FUNCTION 5 .a1 r B a 1 0 0 0 0 6 0 0 0 6 0 0 0 4 0 0 0 2 0 0 0 O V o V © A 3O0O (lO 8) ZCDO looo DISCHARGE - 5 cf.s J I r Q 8 vuxo 2ooo0 1 6 0 0 0 12000 f-8000 4000 - r e g REACH 1 V (mg/l) Fig 6.3 REACH 2 REACH 3 3 3 7 Greg (mg/l) OBJECTIVE FUNCTION 2 V A 2000 00 5 ; 2000 1000 5 7 3 7 3 Fig 6.4 A DISCHARGE* Z'Jcfs (95%) V 6 6 6 J t BREACH 1 REACH2 REACH 3 2 0 C O O L 9 I 6 0 0 O _ 12000 r s ^ 8000 £ 4 0 0 0 A 0 V A © 2 V I O 3 5 7 3000 (io5; 2000 1000 3 & 7 Creg(mg/L) OBJECTIVE FUNCTION PISCHAPiQE 1 1 37-5efs (39%) K,= 0-3 Vd D = zrn.ya A © A © _ i 1 Fig 6.5 3 5 7 Cre3 (mg/L) - 113 -from a 5% mean monthly summer discharge (5 cfs) to the 95% discharge (25 c f s ) , the t o t a l optimal BOD loading of the r i v e r increased about 100% f o r a C of 5 mg/1 and at the 0.9 compliance l e v e l . This suggests reg a r e l a t i v e l y narrow range within which optimal summer BOD loadings are permissable due to the low d i l u t i o n capacity. At 5cfs, loadings i n the four reaches were 4000 to 5000, 1900 to 2500, 0 to 500 and 0 kg/d res p e c t i v e l y for compliance l e v e l s of 0.8 to 1.0 and a C r e g °f 5 mg/1. This i s compared to 10000 to 11000, 600 to 1000, 0 to 500 and 0 kg/d re s p e c t i v e l y for four reaches at 25 c f s . For reach 2, i t appears that at 5 cfs (5% discharge level) and 25 cfs (95% discharge l e v e l ) r e s p e c t i v e l y the optimal loadings were higher f o r the former. This r e s u l t may be due to a higher l e v e l of v a r i a t i o n i n the optimal loadings generated due to trade-offs between reaches made i n the NLP while seeking the optimal value or due to premature termination of the NLP algorithm at coarse convergence c r i t e r i a . This loadings difference may also r e s u l t from the dynamics of the t i d a l gate at low flows where the gate remains closed for extended periods i n the t i d a l cycle with infrequent openings f o r r i v e r discharge. The closed gate would cause r i v e r waters to back up thus increasing the volume of water i n a l l four reaches. This would provide added d i l u t i o n capacity for BOD loadings with reduced convection a b i l i t y . The r e s u l t s suggest that at low flows the d i l u t i o n e f f e c t s i n reach 2 are greater than the l o s t convection e f f e c t s (which would tend to reduce BOD concentrations) allowing for greater loadings. This phenomena i s reduced at higher flows where convection e f f e c t s dominate. Further i n v e s t i g a t i o n would be required to confirm t h i s speculation. For a l l discharge l e v e l s the optimal BOD loadings decreased i n the upstream d i r e c t i o n as d i l u t i o n c a pacities - 114 -became reduced. Reach 4 is particularly sensitive with negligible optimal loadings. These optimal loading levels would be useful in establishing a BOD permit system whereby agricultural polluters would be allowed prorated (by cultivated area) loadings permits depending on their location along the river. However, such a system is weak at best as an abatement strategy since significant compliance monitoring resources would be required since the pollutant sources are largely non-point and, more importantly, d i f f i c u l t to trace to specific agricultural enterprises. Thus, i t would be d i f f i c u l t to relate actual BOD river concentrations to individual agricultural practices making i t d i f f i c u l t to identify enterprises for on-farm treatment strategies, unless enterprise specific monitoring schemes were established, at higher costs. Beyond such a permit system, other measures would be necessary as part of an overall abatement strategy based on longer-term policies but with short-term policies to accommodate transient shocks to the system. The impact of these various measures (aeration, better land use, flow augmentation, waste treatment) on water quality, as defined by the BOD/DO models developed, could be assessed by using the methodology developed previously, since i t provides a rational framework for determining the optimal loading capacity of the river system which is a pre-requisite to policy formulation. Flow augmentation strategies, based results of figures 6.3 to 6.5, could be designed to redistribute the dilutive capacity of the river or introduce additional capacity from outside the system. Reach 1 could be a possible source of waters for dilution in a form of a feedback dilution system to reaches 2, 3 and 4. The optimal volume of - 115 -diversion would be determined such that the increased BOD l e v e l s i n reach 1, due to upstream loadings, would s t i l l meet the DO regulations. This optimal l e v e l could be determined by modifying the hydrodynamic and BOD/DO models to include transient feedback so as to determine steady-state (where concentrations changes due to additions have s t a b i l i z e d ) DO l e v e l s . An i n t e r - b a s i n t r a n s f e r p o l i c y would be required for s i t u a t i o n s where feedback would generate unacceptably low DO concentrations, based on optimal t r a n s f e r l e v e l s f or addition to, or i n replacement of, i n t r a - b a s i n diversions. 6.3.2 DO Regulatory Limit P r i m a r i l y for reach 1 and the objective function value, a reduction i n C allows for higher BOD loadings. Reaches 2 and 3 reg have a high variance with unclear trends with respect to changes i n C as discussed i n the previous section. BOD loadings would reg increase i n an exponential manner with the reduction of C to 0 reg mg/1 where i n f i n i t e optimal loadings would be ( t h e o r e t i c a l l y ) p o s s i b l e . Within the 3 to 7 mg/1 range, the BOD loadings are approximately l i n e a r with highly uncertain p r e d i c t i v e a b i l i t y outside t h i s range. At 5 cfs the optimal loading for 5 and 7 mg/1 regulatory l i m i t s are within 1000 kg/d i n d i c a t i n g an equivalence of the two i f used as an enforced standard. At higher discharges the s e n s i t i v i t y of the optimal loadings to C i s higher i n d i c a t i n g the need to c o r r e c t l y determine the reg regulatory l i m i t when used for seasonal BOD c o n t r o l . 6.3.3 Compliance Levels The t o t a l BOD loading of the r i v e r increase with decrease i n compliance l e v e l . In e s t a b l i s h i n g a f l e x i b l e DO standard two parameters are important: - 116 -a) the DO regulatory l i m i t ; n Pr(LEVEL>LIMIT)>COMPLIANCE b) the acceptable compliance l e v e l s . ) A DO standard which e x p l i c i t l y regulates compliance l e v e l s would be almost impossible to enforce due to the high resource l e v e l required for f i e l d monitoring. However, the compliance l e v e l can be a useful concept i n the design of standards since assumptions of p o l l u t e r compliance behaviour (conscious or non-conscious) can be tested. For example, a 5 mg/1 standard with an assumed 80% compliance l e v e l would r e s u l t i n an e f f e c t i v e standard of 4 mg/1 (from figure 6.2). Thus, a higher c r e ^ could be set for each month such that the e f f e c t i v e DO l i m i t would be 5 mg/1 at some assumed compliance l e v e l . 6.3.4 Model Parameters The s e n s i t i v i t y of optimal loadings to BOD decay (k^) and lo n g i t u d i n a l dispersion (D) i s shown i n f i g u r e 6.6 and 6.7. As expected, the loadings increase with lower k^ and higher D (see equations 3.17 and 3.13) due to the lower impact on DO concentrations of decaying BOD and lower BOD concentrations due to greater BOD dispersion. However, the BOD loadings were more s e n s i t i v e to changes i n the former parameter compared to the l a t t e r , f o r the same percentage change. This has general implications to data gathering and parameter estimation e f f o r t s , namely that greater resources should be a l l o c a t e d i n the determination of k^ than of D. This i s also l i k e l y the case for determination of k 2 although t h i s was not e x p l i c i t l y tested i n t h i s study. From the v a l i d a t i o n a p p l i c a t i o n to the Delaware Estuary i n chapter 5 the s e n s i t i v i t y of decay parameters to stream temperatures was evident and, though not explored i n t h i s study, would require care f u l f i e l d and laboratory work to e s t a b l i s h appropriate temperature - 117 -Sensitivity of BOD Loadings to Longitudinal Dispersion, BOD Decay and Oxygen Source Levels REACH 1 .3> 5 Q 8 •2JOOCO \2CCO 8000 4000 V A 6 REACH Z REACH 3 6 0 0 0 4000 2X30 A Q 1-0 COMPLIANCE A 0-9 COMPLIANCE V 0&COMPLIANCE -reg ( m q / t ) * 7 3 Fig 6.6 OBJECTIVE FUNCTION 5000 (10*) 2000 ICOO DISCHARGE = 2 5 o f s K i = O-S'/d 15 ft 6 J I 3 3 7 Creg (mfl/0 ?EACHvi woco ^ I2QX) 10000 o r . A 8 ecoo -/toco -A O © R E A C H Z goco 2CC0 5 7 Creg (mg/L) - 9 -R E A C H 3 IS OBJECTIVE FUNCTION 30CO do5; 2CCD ICOO A © 7 3 5 Fig 6.7 DISCHARGE = 25Cfs Ki = 005 >/d 0 = Z f M i % l J I 3 & 7 C n e g & n g / L ) B E A C H 1 REACH2 REACH 3 VST r I § •2O0CO \20CO -A © 0000 -40OO -60D0 •4000 -ZOOD 6 SfcOO 1CCD lOOO JiL OBJECTIVE FUNCTION OI5CttARGE=Z5cf5 A © A © J i 3 5 7 3 C r e g ( m g / L ) 7 3 Fig 6.8 3 5 7 Creg (ing/l) ' - 118 -c o e f f i c i e n t s . This would be of great concern f o r summer low flow conditions where higher temperatures would increase BOD decay rates and thus reduce DO concentrations further. The optimal a l l o c a t i o n of data gathering resources f o r maximum information i s an issue f o r further study. 6.3.5 Oxygen Sources As discussed i n section 6.1.4 a net oxygen source l e v e l of 0.09 mg/1 was assumed f o r the four reaches of the DO model. This was a steady source approximation to diurnal photosynthetic processes and could also be interpreted as incorporating steady a r t f i c i a l aeration. The optimal solutions were tested for an increase of 25% to determine the s e n s i t i v i t y of t h i s assumed source l e v e l and to assess the u t i l i t y of longer-term oxygen addition (compared with ad-hoc emergency aeration). The 25% increase i n net oxygen addition (at 25 cfs) r e s u l t e d i n an average increase of about 10% i n BOD loading i n reach 1 (figure 6.8). This lower response may be due to the low l e v e l of DO addition assumed. Order-of-magnitude changes i n photsynthetic DO addition over the diurnal c y c l e , may be required to represent the great v a r i a t i o n s noted i n f i e l d studies (Moore, 1985a). The incorporation of l a t e summer a l g a l decay into the BOD/DO model would further add to describing t h i s v a r i a t i o n . However, the impact of the average DO loading used i n t h i s study could be used to test designs of steady oxygen del i v e r y systems as opposed to ad-hoc intense aeration with diurnal v a r i a t i o n s . 6.3.6 Other Abatement P o l i c i e s As outlined i n section 6.3.1, a m u l t i - l e v e l abatement strategy i s necessary for broader control of water p o l l u t i o n . The optimal BOD loadings generated i n t h i s study are useful i n comparison with actual - 119 -(though d i f f i c u l t to determine) loadings for e s t a b l i s h i n g and combining the components of an o v e r a l l abatement strategy. Additional to the previous p o l i c i e s o u t l i n e d are other long-term p o l i c i e s are suggested: a) On farm treatment and management: Treatment or reduction at source i s a better p o l i c y r e q u i r i n g a d d i t i o n a l on-farm expense with l i k e l y government subsidy support for an already-supported f i n a n c i a l l y stressed industry. An a l t e r n a t i v e would be to encourage the optimal land a p p l i c a t i o n (temporally and s p a c i a l l y ) of animal or other organic wastes, to minimize t h e i r impact on water q u a l i t y through runoff and groundwater seepage. The optimal loadings generated would provide constraints to non-point seepage models, such as Cappelaere (1978), which simulate the impact of pollutants on groundwater q u a l i t y and r i v e r discharges while modeling the e f f e c t of waste nutrients on crop growth. The establishment of optimal farm c u l t i v a t i o n and a p p l i c a t i o n schedules would be an i n t e r e s t i n g multiobjective optimization problem for further study. b) Land use p o l i c i e s : A g r i c u l t u r a l land use p o l i c i e s for marginal land development based on s o i l q u a l i t y and a v a i l a b i l i t y of i r r i g a t i o n water should i d e a l l y account for water q u a l i t y impacts. The impact of a g r i c u l t u r a l development along reaches 2, 3 and 4 would be reduced by the diversion of r i v e r inflow i n the form of land runoff, d i t c h flows and seepage towards reach 1 with i t s greater d i l u t i o n capacity and/or the use of instream aeration to compensate for oxygen demands i n sub-reaches during c r i t i c a l low DO periods such as a l g a l d i e - o f f and spawning. - 120 -Throughout t h i s study, various model parameters (such as k^, k 2/ D) and inputs (such as discharge, temperature) were assumed to be dete r m i n i s t i c . Ward and L o f t i s (1983) stressed the balance between deterministic and stochastic view of water q u a l i t y modeling and management. This stochastic nature could have been incorporated into the models developed i n t h i s study, based on simple moments of normal and log-normal d i s t r i b u t i o n s , by random generation of model parameter values possible using Markov techiques. This would have embelished the joint-chance constraint NLP formulation of section 3.1 and would be an i n t e r e s t i n g enhancement for further study. The discussions above were based on r e s u l t s of t h i s study derived from an unvalidated BOD/DO model. Even with s u f f i c i e n t v a l i d a t i o n the a p p l i c a b i l i t y of the model would s t i l l be tenuous since p o l i c i e s bases on l i m i t e d dimensions of a system do not contain the r e q u i s i t e v a r i e t y to m u l t i - v a r i a b l y control the complex ecosystem processes of organic p o l l u t i o n i n dynamic r i v e r s . - 121 -7. GENERAL DISCUSSION 7.1 Hydrodynamic Model A good model of the hydrodynamics i s c r u c i a l to the accurate simulation of the water q u a l i t y i n a dynamic r i v e r . The model should account for the major sources and sinks of water i n the r i v e r system including t r i b u t a r y flow, groundwater flow ( i f s i g n i f i c a n t ) , i r r i g a t i o n (uptake and return flow) and runoff. The model used i n t h i s study did not simulate these l a t t e r processes but accounted only f or upstream inflows and downstream t i d e s , which appears to adequately match h i s t o r i c a l records. There are more e f f i c i e n t f i n i t e - d i f f e r e n c e schemes and algorithms for the sol u t i o n of the equations of continuity and momentum (Mauersberger, 1983), allowing for shorter model execution times, which can be s i g n i f i c a n t on a microcomputer. The extension of the model i n t o two dimensions may be of research i n t e r e s t , e s p e c i a l l y f o r the p r e d i c t i o n of po l l u t a n t plumes. Another extension of great i n t e r e s t would be the modeling of the e s t u a r i a l conditions of Mud Bay with multiple pollutant loadings from both the Nicomekl and Serpentine r i v e r s . This would help assess the d i l u t i o n and mixing capacities of the bay. The use of stochastic hydrographs as upstream boundary condition generators (as done by Biswas and Reynolds, 1973) combined with computed tides (from harmonic equations) could be useful i n p r e d i c t i n g seasonal or annual stream flow and hence, water q u a l i t y v a r i a t i o n s . This would of value f or longer-term decision horizons requiring, however, greater computing power. This extension would also - 122 -f a c i l i t a t e Monte Carlo simulation of various stochastic model parameters. In reviewing the hydrodynmics of the Nicomekl River during summer low flow conditions, i t appears that a steady-state flow model may be applicable without a great loss of information and reduce computing e f f o r t s . This i s of use i f the BOD loading l i m i t s were to be regulated based on an i n f l e x i b l e DO standard based on low flow (high temperature) extremes. A f l e x i b l e regulatory system would necessarily have to accommodate dynamic r i v e r processes. 7.2 BOD/DO Model The BOD/DO model used i n t h i s study was based on a compartmentalization of the r i v e r f o r the purposes of pol l u t a n t mass balances. More e f f i c i e n t schemes and algorithms also e x i s t for such models with extensions into broader water q u a l i t y / e c o l o g i c a l formulations which provide greater number of dimensions for regulatory control and p o l i c y evaluation (Mauersberger, 1983). The a p p l i c a t i o n of estimation techniques could be used to amplify the l i m i t e d f i e l d data a v a i l a b l e i n i d e n t i f y i n g possible stochastic, steady-state BOD/DO models (Damaskos and Papadopoulos, 1983) with p r o b a b i l i s t i c v a r i a b l e s such as stream flow (to r e f l e c t the t i d a l and hydrograph dynamics); stream temperature (and hence model" parameters: k^, k 2/ D); i n i t i a l and boundary conditions (for dynamic inputs); and measurement (for BOD and DO t e s t i n g , the former being o f f - l i n e and more uncertain); (Beck, 1981,1983,1985). For the Serpentine r i v e r , t h i s approach would have been i n t e r e s t i n g , not only i n model i d e n t i f i c a t i o n but also i n terms of real-time DO c o n t r o l , as i s currently underway (Mavinic, 1986). - 123 -A simpler steady-state model would also have been easier to incorporate into the NLP formulation and, again, would require less time for model execution on the microcomputer. The r e s u l t s would have to be more c a r e f u l l y interpreted, e s p e c i a l l y for short-term a p p l i c a t i o n s . 7.3 Uncertainty of Measured BOD As mentioned i n section 2.2.1, the use of BOD as a measure of water q u a l i t y has an inherent measurement problem. Standard t e s t s measure BOD as a 5-day BOD or the l e v e l of BOD a f t e r a 5-day incubation. The BOD used i n the models formulated are based on the ultimate BOD l e v e l , or BOD , . which i s often modeled as u l t BOD^ = (1 - e" K L n)BOD _ n u l t where i s a laboratory determined decay parameter. Thus, the optimal BOD l e v e l s discussed further i n t h i s chapter have to be converted to BOD^ so as to be used for actual f i e l d control purposes. P h e i f f e r et a l (1976) found that laboratory errors i n the measurement of BOD made i t a non-ideal water q u a l i t y parameter without s u f f i c i e n t t e s t i n g for s t a t i s t i c a l v a l i d i t y . This v a r i a t i o n could be modeled as measurement noise using estimation techniques as outlined i n section 2.2.4. A l t e r n a t i v e l y , more consistent real-time water q u a l i t y i n d i c a t o r s such as t o t a l organic carbon (TOC) could be used to estimate organic oxygen demand. 7.4 Data Gathering Data gathering i s an expensive exercise requiring optimization also so as to obtain the maximum amount of information subject to some f i e l d budget. The frequency, i n t e n s i t y and q u a l i t y of information gathered i s d i f f e r e n t for mathematical model b u i l d i n g purposes and f o r real-time forecasting and c o n t r o l , with greater e f f o r t required for the - 124 -former. An i n t e r e s t i n g methodology for network design based on modern concepts of communication theory was developed by Husain (1980) applicable to watersheds of the si z e of the Nicomekl/Serpentine basin. The present data gathering e f f o r t s need to be r a t i o n a l i z e d on the basis of experimental designs for mathematical model b u i l d i n g and v a l i d a t i o n , with f i e l d experiments rel a t e d s p e c i f i c a l l y to convection, dispersion and decay processes of water p o l l u t i o n for more water q u a l i t y measures than BOD. 7.5 NLP Formulation There are possible s i m p l i f i c a t i o n s of the NLP formulation i n section 3.1. These are: i ) the decoupling of the j o i n t constraints for derivation of deterministic equivalents while s t i l l incorporating the dynamic DO model. This would then r e f l e c t separate allowable r i s k l e v e l s for the d i f f e r e n t reaches considered incorporating r e a c h - s p e c i f i c information. i i ) the NLP formulation could be approximated by an i t e r a t i v e LP method where l i n e a r constraints could be derived from the regression of constraint contours i n appropriate constraint dimensions (such as reach 2 and reach 3 i n section 6.2). The objective function would be evaluated i n various d i s c r e t e nonlinear dimensions and the maximum obtained. This approach, however, would require g r i d searches i n i d e n t i f y i n g l i n e a r dimensions which may be p r o h i b i t i v e f or higher dimensioned problems. 7.6 NLP Algorithm The Hooke-Jeeves algorithm was adequate for the optimization problem formulated, though a fa s t e r converging algorithm (perhaps Powell's) may have have marginally reduced the elapsed time for ) - 125 -searching for the optimal so l u t i o n . The computations involved i n the DO model were by far the l i m i t i n g factor i n running the NLP algorithm, so that e f f o r t s i n streamlining the model scheme would be more productive. A d d i t i o n a l l y , a comparison of the s t a b i l i t y of t h i s scheme with other d i r e c t search techniques (such as Powell, Nelder-Mead and Rosenbruck) would also have been of great value i n i d e n t i f y i n g the better algorithm for solving the problem formulated i n section 3.1. 7.7 Computer Resource The COMPAQ 286 microcomputer provides l i m i t e d resources for the implementation of the methodology proposed i n section 3.3.1 at no cost to research or hard-dollar mainframe accounts. With the previously mentioned modifications to increase model e f f i c i e n c y , the execution times for a four-reach, seven-day problem may be reduced s i g n i f i c a n t l y . Larger problems, less than 640 KB i n s i z e , would take longer to run providing greater j u s t i f i c a t i o n for mainframe or 32-bit microcomputer use. The microcomputer i s a good research and model development t o o l for small- to medium-sized s t r a t e g i c analyses (monthly, seasonal, annual) of water q u a l i t y management problems as judged from the l i m i t e d a p p l i c a t i o n to non-linear optimization i n t h i s study. However, due to i t s l i m i t e d processing c a p a b i l i t i e s , the machine may not be appropriate for real-time control or analysis where short (diurnal) p r e d i c t i o n or optimization may be important for timely decision making. 7.8 U t i l i t y of Proposed Methodology to Management P r a c t i s e The methodology presented in t h i s study with enhancements (model v a l i d a t i o n , use of robust NLP algorithms, use of a steady-state pollutant model) could be used by management to r a t i o n a l i z e solutions - 126 -for f i e l d problems and formulating short- and long-term abatement p o l i c i e s . The (near) optimal solutions generated by the methodology provide guideline BOD loading l e v e l s which i n comparison with actual loadings would help i d e n t i f y the nature of p o l l u t i o n control technique required. - 127 -8. CONCLUSIONS 1) Using a one-dimensional hydrodynamic model i t was possible to c o r r e c t l y p r e d i c t the stage and discharge at various cross-sections of a t i d a l r i v e r . 2) Using a one-dimensional BOD/DO model coupled with 1) above i t was possible to describe and, with less confidence due to the lack of s i t e s p e c i f i c data for v a l i d a t i o n , predict the r e s u l t i n g concentrations due to the processes of convection, dispersion and decay. 3) Using a non-linear programming algorithm coupled with the DO model 2), i t was possible to maximize the t o t a l BOD loadings of a t i d a l r i v e r with various DO regulatory l i m i t s and r i s k l e v e l s of v i o l a t i o n . 4) A powerful microcomputer i s an e f f e c t i v e t o o l f o r one-dimensional hydrodynamic and water q u a l i t y modeling. However, execution times may become p r o h i b i t i v e compared with mainframe computers but at zero cost. The methodology proposed i n t h i s study i s s u i t a b l e f or 'smaller-sized' problems where the l i m i t e d computing resources are not a serious constraint. 5) The methodology presented appears to be a v i a b l e method of determining the optimal pollutant loadings of a t i d a l r i v e r f o r s t r a t e g i c decision-making. However extensive model i d e n t i f i c a t i o n should be conducted to j u s t i f y the choice of water q u a l i t y models. 6) The water q u a l i t y data gathering on the Nicomekl/Serpentine r i v e r s should be based on a more complete experimental design to allow for - 128 -water q u a l i t y model development and v a l i d a t i o n , and for longer-term po l l u t a n t and DO con t r o l requirements. - 129 -REFERENCES Ages, A. (1973). A Numerical Model of V i c t o r i a Harbour to P r e d i c t T i d a l  Responses to Proposed Hydraulic Structures. P a c i f i c Marine Science Report No. 73-3, Marine Sciences Branch, Environemnt Canada. Aiba, S and H. Ohtake (1977). Simulation of P04-P Balance i n a Shallow  and P o l l u t e d River. Water Research, Vol 11, 159-164. Ambrose, R (1986). Personal Communication. USEPA, Athens, Georgia. Arbabi, M., J. Elzinga and C S . Revelle (1974). The Oxygen Sag Equation:  New Properties and a Linear Equation f o r the C r i t i c a l D e f i c i t . Water Resour. Res., 10, 921-929. Bansal, M.K. (1976). Dispersion Model f o r an Instantaneous Source of P o l l u t i o n i n Natural Streams and i t s A p p l i c a b i l i t y to the B i g Blue  River (Nebraska). In W.R. Ott 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 335-339. Barr, L.J. and D.E. Rekstein (1984). Nicomekl-Serpentine Water Supply Study. Memorandum, Water Management Management Branch, Ministry of Environment, Province of B r i t i s h Columbia. Bayer, M.B. (1972). A Nonlinear Mathematical Programming Model f o r Water  Quali t y Management. In A.K. Biswas (ed) 'International Symposium on Modeling Techniques i n Water Resources Systems', Environment Canada, Ottawa, 341-351. Beck, M.B. (1981). Operational Water Quality Mangagement: Beyond Planning and Design. International I n s t i t u t e for Applied System Analysis, Executive Report 7. Beck, M.B. (1983). A Procedure for Modeling and S e n s i t i v i t y A n a l y s i s ,  C a l i b r a t i o n and V a l i d a t i o n . In G.T. Orlob (ed) "Mathematical Modeling of Water Quality: Streams, Lakes, and Reservoirs'. I n s t i t u t e f o r Applied System Analysis and John Wiley and Sons. Beck, M.B. (1985). Water Quality Management: A Review of the Development  and A p p l i c a t i o n of Mathematical Models. I n s t i t u t e for Applied System Analysis, Laxenburg. Beckers, CV. , P.E. Parker, R.N. Marshall and S.G. Chamberlain (1976). RECEIV-II: A Generalized Dynamic Planning Model f o r Water Qual i t y  Management. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 344-349. Beer, S (1966). Decision and Control. John Wiley, Chichester. B e l l a , D.A. (1970). Dissolved Oxygen Va r i a t i o n s i n S t r a t i f i e d Lakes. J. San. Eng. Div., ASCE, SA5, 1129-1146. - 130 -B e l l a , D.A. and W.E. Dobbins (1967). F i n i t e - D i f f e r e n c e Modeling of River  and Estuary P o l l u t i o n . Proceedings of the National Symposium on Estuarine P o l l u t i o n , Stanford University, 612-645. B e l l a , D.A. and W.E. Dobbins (1968). F i n i t e - D i f f e r e n c e Modeling of Stream P o l l u t i o n . J. San. Eng. Div., ASCE, v o l . 94, SA5, 995-1016. B e l l a , D.A. and W.J. Grenney (1970). F i n i t e - D i f f e r e n c e Convection  Errors. J. San. Eng. Div. ASCE, SA6, 1361-1374. B e l l , K.M. and A.S. Papadopoulos (1979). A Stochastic Model f o r BOD and  DO when the Discharged Pollutants and the Flow of the Stream are  Random Quantities. Intern. J. Environ. Studies, 11, 37-42. Bishop, A.B. and W.J. Grenney (1976). Coupled Optimization-Simulation  Water Qual i t y Model. J. Env. Eng Div. ASCE, v o l . 102, EE5, 1071-1086. Biswas, A.K. and P.J. Reynolds (1973). Mathematical Models f o r the  Planning and Management of the Saint John River Systems. International Symposium on River Mechanics, IAHR, Bangkok, Thailand, 415-428. Bourque, S. and G. Herbert (1982). A Preliminary Assessment of Water Quality and B i o t a i n the Serpentine and Nicomekl Rivers and Mahood  Creek. Reg. Prog. Rep. 82-17, Environmental Protection Service, Dept. of Environment. Bracken, J. and G.P. McCormick (1968). Deterministic Non-Linear Equivalents f o r Stochastic Programming Problems. In 'Selected Applications of Nonlinear Programming', Wiley. Burges, S.J. and D.P. Lettenmaier (1975). P r o b a b i l i s t i c Methods i n  Stream Qua l i t y Management. Water Res. B u l l . , Vol. 11, No. 1, 115-130. Burn, D.H. and E.A. McBean (1985). Optimization Modeling of Water Quality i n an Uncertain Envoronment. Water Resour. Res., 21(7), 934-940. Chang, S. and W.W-G. Yeh (1973). Optimal A l l o c a t i o n of A r t i f i c i a l Aeration along a P o l l u t e d Stream Using Dynamic Programming. Water Resour. Res, 9(5), 985-997. Chichester, F.D. (1976). Ap p l i c a t i o n of M u l t i l e v e l Control Techniques to  Classes of D i s t r i b u t e d Parameter Plants. Unpublished D.Eng. S c i . Thesis, Dept. of E l e c t r i c a l Engineering, New Jersey I n s t i t u t e of Technology, Newark, New Jersey, Chow, V.T. (1959). Open Channel Flow. McGraw-Hill Book Co. Inc. - 131 -Covar, A.P. (1976). S e l e c t i n g the Proper Reaeration C o e f f i c i e n t f o r Use i n Water Qual i t y Models. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n n a t t i , 340-343. Crean, P.B. (1983). The Development of Rotating, Non-Linear Numerical  Models (GF2, GF3) Simulating Barotropic Mixed Tides i n a Complex  Coastal System Located between Vancouver Island and the Mainland. Can. Tech. Rep. Hydrogr. Ocean S c i : 31, 65p. Currie, M.V. (1984). I r r i g a t i o n Withdrawals i n the Fraser Delta Planning  Unit. Water Management Branch, Ministry of Environment, Province of B r i t i s h Columbia. Damaskos, S.D. and A.S Papadopoulos (1983). A General Stochastic Model  fo r P r e d i c t i n g BOD and DO i n Streams. J. Environ. Management, 15, 229-237. Davidson, B. and R.W Bradshaw (1970). A Steady State Optimal Design of  A r t i f i c i a l Induced Aeration i n Polluted Streams by the Use of  Pontryagin's Minimum P r i n c i p l e . Water Resour. Res., 6(2), 383-397. Davis, J.C. (1975). Minimal Dissolved Oxygen Requirements of Aquatic  L i f e with Emphasis on Canadian Species: A Review. J. F i s h . Res. Board Can., Vol 32, No. 12, 2295-2332. Deininger, R.A. ( 1972). Minimum Cost Regional P o l l u t i o n Control System. In A.K. Biswas (ed) 'International Symposium on Mathematical Modeling i n Water Resources System', Environment Canada, Ottawa, 352-361 Delos, C.G. (1976). Mathematical Model of a Great Lakes Estuary. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 115-119. Dobbins, W.E. (1964). BOD and Oxygen Relationships i n Streams. J . San. Eng. Div., ASCE, Vol. 90, SA3, 53-78. Dresnack, R. and W.E. Dobbins (1968). Numerical Analysis of BOD and DO  P r o f i l e s . J. San. Eng. Div., ASCE, Vol. 94, SA5, 789-807. Dronkers, J . J . (1969). T i d a l Computations f o r Rivers, Coastal Areas and  Seas. J. Hyd. Div., ASCE, Vol. 95, HY1, 29-77. Duffy, R.G. and A.B. Clymer (1976). Water Modeling i n Ohio EPA. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 517-521. Dunne, T. and L.B. Leopold (1978). Water i n Environmental Planning. W.H. Freeman and Company, San Francisco. Environment Canada (1985). H i s t o r i c a l Streamflow Summary: B r i t i s h  Columbia to 1984. Inland Waters Directorate, Water Resources Branch, Ottawa. - 132 -Ford, D.E., K.W. Thornton, A.S. Lessem and J.L. Norton (1980). A Water  Quality Management Model f o r Reservoirs. Proceedings of the Symposium on Surface Water Impoundments, ASCE, June 2-5, Minneapolis, Minnesota, 624-633. Foreman, M.G.G. (1977). Manual for T i d a l Heights Analysis and P r e d i c t i o n . P a c i f i c Marine Science Report 77-10, I n s t i t u t e of Ocean Sciences, P a t r i c i a Bay, V i c t o r i a , B.C. 101pp. Friedman, R., C. A n s e l l , S. Diamond and Y.Y. Haimes (1984). The Use of  Models f o r Water Resources Management, Planning and P o l i c y . Water Resour. Res. 20(7), 793-802. Gorelick, S.M. (1982). A Model for Managing Sources of Groundwater  P o l l u t i o n . Water Resour. Res., 18(4), 773-781. Gorelick, S.M. and I. Remson (1982). Optimal Dynamic Management of Groundwater Pol l u t a n t Sources. Water Resour. Res., 18(1), 71-76. Gourishankar, V. and R.L. Lawson (1975). Optimal Control of Water P o l l u t i o n i n a River Stream. Int. J. Systems S c i . , 6(3), 201-216. Gourishankar, V. and K. Ramar (1977). Design of a Water Q u a l i t y C o n t r o l l e r by Pole Assignment. IEEE Transactions on Systems, Man, and Cybernetics, Vol. SMC-7, No. 2, 109-111. Gourishankar, V. and M.A. Lawal (1978). A D i g i t a l Water Q u a l i t y C o n t r o l l e r f o r P o l l u t e d Streams. Int. J. Systems S c i . , 9(8), 899-919. Gromiec, M.J., D.P. Loucks and G.T. Orlob (1983). Stream Q u a l i t y  Modeling. In G.T. Orlob (ed) 'Mathematical Modeling of Water Quality: Streams, Lakes, and Reservoirs', IIASA and John Wiley and Sons. Haimes, Y.Y. (1976). H i e r a r c h i c a l Analysis of Water Resources Systems, McGraw-Hill Book Co. Inc, Himmelblau, D.M. (1972). A Uniform Evaluation of Unconstrained Optimization Techniques. In F.A. Lootsma (ed) 'Numerical Methods fo r Non-Linear Optimization', Academic Press, London, 69-77. Hung, H.D., A. Hossain and T.P. Chang (1976). Stream Modeling and Waste  Load A l l o c a t i o n . In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 129-132. Husain, T. (1979). Shannon's Information Theory and Hydrologic Network  Design and Estimation. Unpublished Ph.D. Thesis, Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver. - 133 -Ikeda, S. and Y. Sawaragi (1980). Nonlinear Forecasting Models of Water Flow and Q u a l i t y by H e u r i s t i c Self-Organization Methods. In E.F. Wood (ed) 'Real-Time Forecasting/Control of Water Resource Systems', Pergamraon Press. Jacoby, H.D. and D.P. Loucks (1972). Combined Use of Optimization and  Simulation Models i n River Basin Planning. Water Resour. Res., 8(6), 1401-1414. J a i n , R.K. and M.M. Denn (1976). Short-Term Regulation of BOD Upsets i n  an Estuary. J . Dynamic Systems, Measurement and Control, Trans, of the ASME, March, 30-31. Jamal, I.B. (1980). Stochastic River Modeling to Estimate Flood P r o b a b i l i t i e s of a T i d a l River. Unpublished M.A.Sc. Thesis, Dept. of C i v i l Engineering, University of B r i t i s h Columbia. 78pp. Johnson, P.A., M.W. Lorenzen and W.W. Waddel (1976). A Multi-Parameter  Estuary Model. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 111-114. Jorgensen, S.E. (1983). Modeling the E c o l o g i c a l Process. In G.T Orlob (ed) 'Mathematical Modeling of Water Quality: Streams, Lakes, and Reservoirs', IIASA and John Wiley and Sons. Kolbin, V.V. (1976). Stochastic Programming. Theory and Decision Lib r a r y , V ol. 14, D. Reidal Publishing Company. 195p. Kuester, J.L. and J.H. Mize (1973). Optimization Techniques with  FORTRAN. McGraw-Hill Book Company. Lasdon, L.S. (1970). Optimization Theory f o r Large Systems. The Macmillan Company, New York. Lawal, M.A. (1977). Design of Discrete-Time Water Q u a l i t y C o n t r o l l e r s  f o r a P o l l u t e d River System. Unpublished M.Sc. Thesis, Dept. of E l e c t r i c a l Engineering, University of Alberta, Alberta. Lohani, B.N. and N.C Thanh (1978). Stochastic Programming Model f o r Water Qual i t y Management i n a River. J. Water P o l l u t . Contr. Fed., 50(9), 2175-2182. Loucks, D.P. (1983). Models f o r Management A p p l i c a t i o n . In G.T Orlob (ed) 'Mathematical Modeling of Water Quality: Streams, Lakes, and Reservoirs', IIASA and John Wiley and Sons. Mauersberger, P. (1983). General P r i n c i p l e s of Deterministic Water  Quali t y Modeling. In G.T. Orlob (ed) 'Mathematical Modeling of Water Quality: Streams, Lakes, and Reservoirs', IIASA and John Wiley and Sons. Mavinic, D. (1986). Personal Communication. Department of C i v i l Engineering, University of B r i t i s h Columbia. - 134 -M i l l e r , B.L. and H.M. Wagner (1965). Chance Constrained Programming with  J o i n t Constraints. Operation Res., 13, 930-945. Ministry of Environment, Province of B.C. (1984). Water Management  Information System: P i l o t Project. Water Management Branch. Moore, B. (1985a). Serpentine River Dissolved Oxygen Survey 1981-83. Environmental Section, Waste Management Branch, Mini s t r y of Environment, Province of B r i t i s h Columbia. Moore, B. (1985b). Personal Communication. Waste Management Branch, Ministry of Environment, Province of B r i t i s h Columbia. Olgac, N.M., CA. Cooper and R.W. Longman (1976). Optimal A l l o c a t i o n of  Measurement and c o n t r o l Resources with A p p l i c a t i o n to River  Depollution. IEEE Transactions on Systems, Man and Cybernetics, Vol SMC-6, No 4, 377-384. Ozgoren, M.K. , R.W. Longman and CA. Cooper (1975). Optimal Control of a  Di s t r i b u t e d Parameter Time-Lagged River Aeration System. J. of Dynamic Systems, Measurement, and Control, Trans, of ASME, Vol. 97, Series G, No. 2, June, 164-171. Padgett, W.T. and S.D. Durham (1976). A Random D i f f e r e n t i a l Equation A r i s i n g i n Stream P o l l u t i o n . J . Dynamic Systems, Measurement, and Control, Trans, of ASME, Vol. 98, 32-36. Padgett, W.T. and A.S. Papadopoulos (1979). Stochastic Models f o r BOD  and DO i n Streams. Ecolo g i c a l Modeling, 6, 289-303. Padgett, W.T., G. Schultz and C P . Tsokos (1977). A Stochastic Model f o r  BOD and DO i n Streams when Pollutants are Discharged over a  Continuous Stretch. Intern. J. Environ. Studies, 11, 45-55. Paul, J. (1986). Personal Communication. USEPA, Naragansa, Rhode Island. P h e i f f e r , T.H., L.J. Clark and N.L. Lovelace (1976). Patuxent River Basin Model Rates Study. In W.R Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 133-138. Putz, G. , R. Gerard and D.W. Smith (1983). E f f l u e n t D i l u t i o n and Decay,  Slave River Downstream of Fort Smith, N.W.T. Environmental Engineering Technical Report Series 83-1, Dept. of C i v i l Engineering, University of Alberta, Alberta. Ramm, A.E. (1976). A Mathematical Model of Dissolved Oxygen i n the Lower  Cuyahoga River. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 101-105. Rickert, D.A., W.G. Hines and S.W. McKenzie (1976). Planning Implication of Dissolved Oxygen Depletion i n the Willamette River,  Oregon. In W.R Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation", C i n c i n a t t i , 62-68. 135 -Ri n a l d i , S., R. Soncini-Sessa, H. Stehfest and H. Tamura (1979). Modeling and Control of River Q u a l i t y . McGraw-Hill Inc., London, 380pp. Sageev, A. (1986). Slug Test Analysis. Water Resour. Res., 22(8), 1323-1333. Sargent. D.H. (1976). Waste A l l o c a t i o n i n the B u f f a l o (New York) River  Basin. In W. R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 126-128. Sch o l l , J.E. and R.L. Wycoff (1981). Continuous DO Simulation at S p r i n g f i e l d Missouri. J. Env. Eng., ASCE, Vol.107, EE1, 69-82. Schultz International Ltd. (1983). Development of a Reservoir Water Quality Model f o r Reservoirs i n B r i t i s h Columbia. A report prepared for B.C. Hydro and Power Authority, Province of B r i t i s h Columbia, Project B114.3.39. Sengupta, J.K. (1972). Stochastic Programming. North-Holland Publishing Company, Amsterdam. Shastry, J. S. , L.T. Fan and L.E. Erickson (1973). Nonlinear Parameter  Estimation i n Water Qual i t y Modeling. J. Env. Eng. Div., ASCE, V o l . 99, EE3, 315-331. Singh, M.A. (1977). Dynamical H i e r a r c h i c a l C o n t r o l . North-Holland Publishing Company, 257pp. Spear, R.C. and G.M. Hornberger (1983). Control of DO l e v e l i n a River  Under Uncertainty. Water Resour. Res., 19(5), 1266-1270. Tabak, D. and B.C. Kuo (1970). Optimal Control by Mathematical  Programming. P r e n t i c e - H a l l , New Jersey. Tamura, H (1974). A Discrete Dynamic Model with D i s t r i b u t e d Transport  Delays and i t s H i e r a r c h i c a l Optimization f o r Preserving Stream  Quality. IEEE Trans, on Systems, Man, and Cybernetics, Vol. SMC-4, No. 5, September, 424-431. Tamura, H. and T. Kawaguchi (1980). Real-Time Estimation of Orders and  Parameters of Distibuted-Lag Models of River Q u a l i t y . In E.F. Wood (ed) 'Real-Time Forecasting/Control of Water Resource Systems', Pergamraon Press. Tarrasov, V.J., H.J. P e r l i s and B. Davidson (1969). Optimization of a  Class of River Aeration Problems by the use of M u l t i v a r i a b l e D i s t r i b u t e d Parameter Control Theory. Water Resour. Res., 5(3), 563-573. Taylor, A.C. and R.W Judy (1972). Waste Abatement Control by Rent Allotment. In A.K. Biswas 'International Symposium on Mathematical Modeling Techniques i n Water Resources Systems', Environment Canada, Ottawa. 136 -Thayer, R.P. and R.G. Krutchkoff (1967). Stochastic Model f o r BOD and DO  i n Streams. J. San. Eng. Div.,ASCE, Vol. 93, SA3, 59-72. Thompstone, R.M., D.K. Sen and R. Divi (1979). Hydrologic and Optimization Techniques for Operating Multi-Reservoir Water  Resources Systems. Fourth National Hydrotechnical Conference on River Basin Management, CSCE, Vancouver, B.C., 132-150. Tubbs, L.J. and D.A. Haith (1981). Simulation Model f o r A g r i c u l t u r a l  Non-Point-Source P o l l u t i o n . JWPCF, 53(9), 1425-1433. Tucci, C.E.M. and Y.H Chen (1981). Unsteady Water Q u a l i t y Model f o r River Network. J. Water Res. Plan, and Man., ASCE, Vol.107, WR2, 477-493 Van der Gulik, T. (1983). Fraser River S t r a t e g i c Plan. Memorandum, A g r i c u l t u r a l Engineering Branch, Ministry of A g r i c u l t u r e and Food, Province of B r i t i s h Columbia. Van der Gulik, T. (1986). Personal Communication. Mini s t r y of Ag r i c u l t u r e and Food, Province of B r i t i s h Columbia. Walsh, G.R. (1979). Methods of Optimization. John Wiley and Sons, London Ward, R.C. and J.C L o f t i s (1983). Incorporating the Stochastic Nature of  Water Qual i t y i n t o Management. JWPCF, 55(4), 408-414. Weatherbe, D.G. (1976). A Dynamic Water Qual i t y Simulation Model f o r the  Thames River. In W.R Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 330-334. Wen, C-G., J-F. Kao, L.K. Wang and M.H. Wang (1982). Determination of  S e n s i t i v i t y of Water Qual i t y Parameters f o r Stream P o l l u t i o n  Control. J. of Environ. Management, Vol. 14, 17-34. Wood, E.F. (ed) (1980). Real-Time Forecasting/Control of Water Resource  Systems. Selected papers from an IIASA Workshop, 1976. Pergammon Press. Yakowitz, S. (1982). Dynamic Programming App l i c a t i o n s i n Water  Resources. Water Resour. Res., 18(4), 673-696. Yeasley, J.R (1976). The A p p l i c a t i o n of a Steady-State Water Qual i t y  Model to the Permit Writing Process, Lake Milner, Idaho. In W.R. Ott (ed) 'Proceedings of the EPA Conference on Environmental Modeling and Simulation', C i n c i n a t t i , 780-783. Yeh, W.W-G. (1985). Reservoir Management and Operations Models; A State-of-the-Art Review. Water Resour. Res., 21(12), 1797-1818. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0096653/manifest

Comment

Related Items