UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Water quality modelling in estuaries Joy, Christopher Stewart 1974

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
UBC_1974_A1 J69.pdf [ 8.79MB ]
[if-you-see-this-DO-NOT-CLICK]
Metadata
JSON: 1.0062998.json
JSON-LD: 1.0062998+ld.json
RDF/XML (Pretty): 1.0062998.xml
RDF/JSON: 1.0062998+rdf.json
Turtle: 1.0062998+rdf-turtle.txt
N-Triples: 1.0062998+rdf-ntriples.txt
Original Record: 1.0062998 +original-record.json
Full Text
1.0062998.txt
Citation
1.0062998.ris

Full Text

WATER QUALITY MODELLING IN ESTUARIES by CHRISTOPHER STEWART JOY B.Eng., Mohash University, 1968 M.App.Sci., Monash University, 1972 A THESIS. SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF•PHILOSOPHY in the Department of Civil Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA January, 1974 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Q U\^kSf (jrS&^oJic,y> The University of British Columbia Vancouver 8. Canada Date ABSTRACT A tidally varying and a tidally averaged mass transport model are applied to the Fraser River Estuary to investigate the significance of tidal effects on the concentrations resulting from assumed effluent discharges. The tidally averaged model is due to Thomann [1963]. The tidally varying model is developed from first principles. A hydrodynamic model was used to determine the tidally induced temporal variation in the longitudinal velocity and cross-sectional area along the estuary. All models are "mathematical" and one-dimensional. Finite difference techniques are used to solve the underlying partial differential equations of all three models. The problems of stabil ity and numerical dispersion are examined. Numerical dispersion is seen to result from the solution of the mass transport equation over a fixed space grid rather than along the advective characteristics. Advantages of solving the equation along the characteristics are: no numerical dispersion; the advective and dispersive transport processes are usefully separated; lateral dispersion can be partially assessed with a one-dimensional model; and time dependent behaviour in coefficient of longitudinal dispersion can be taken into account. The tidally varying flows along the estuary are seen to cause a variation in the initial dilution of a discharged effluent. This, together with the effects of tidal flow reversal produces spikes in the concentration profile along the estuary. The concentration of these spikes is then re duced by the dispersion process, the peak concentration during the first two ii iii tidal cycles being sensitive to the form and magnitude of the coefficient of longitudinal dispersion. Time dependent variations in this coeffi cient are considered. The effect of the lateral dispersion process on the estimated concentrations is also considered and secondary flows are tentatively explained in terms of the generation and advection of vor-ticity. The predicted peak tidally varying concentration was found to be significantly greater than the tidally averaged value. TABLE OF CONTENTS Page LIST OF TABLES vii LIST OF FIGURES viiLIST OF SYMBOLS xCHAPTER INTRODUCTION. ... 1 1. PRELIMINARY CONSIDERATIONS 5 1.1 THE ONE-DIMENSIONAL MASS TRANSPORT EQUATION 5 1.2 DETERMINATION OF PARAMETERS 7 1.3 THE EFFECTS OF THE TIDE ON MASS TRANSPORT PROCESSES 9 1.4 ACCURACY OF A ONE-DIMENSIONAL MODEL 12 2. LITERATURE REVIEW 17 2.1 ANALYTICAL SOLUTIONS 12.2 NUMERICAL SOLUTIONS 21 2.2.1 Finite Difference Solutions 22.2.2 "Box Model" Solutions 24 2.2.3 Finite Element Solutions 5 2.3 PHYSICAL AND ANALOGUE MODEL SOLUTIONS 27 2.3.1 Physical Model Solutions 22.3.2 Analogue Model Solutions 8 2.4 STOCHASTIC SOLUTIONS 29 2.5 SUMMARY 32 3. A DESCRIPTION OF THE HYDRODYNAMIC AND MASS TRANSPORT MODELS ... 34 3.1 THE HYDRODYNAMIC MODEL 33.1.1 The Hydrodynamic Equations 333.1.2 Assumed Quasi-Steady Hydraulic Conditions 38 3.1.3 The Model Estuary of the Hydrodynamic Equations ... 39 3.2 THE TIDALLY VARYING MASS TRANSPORT MODEL 33.2.1 Method of Solution 33.2.2 The Model Estuary 42 3.3 TIDALLY AVERAGED MASS TRANSPORT MODEL 44 3.3.1 Method of Solution3.3.2 The Model Estuary 7 iv V CHAPTER Page 4. VERIFICATION OF THE HYDRODYNAMIC AND MASS TRANSPORT MODELS. ... 48 4.1 THE HYDRODYNAMIC MODEL 44.1.1 Data Available4.1.2 Water Surface Elevations for Low Flows ........ 49 4.1.3 Water Surface Elevations for High Flows 44.1.4 Verification for the Conditions of January 24, 1952 . 54 4.2 THE TIDALLY VARYING MASS TRANSPORT MODEL. 65 4.2.1 Advective Transport 64.2.2 Dispersive Transport 7 4.3 THE TIDALLY AVERAGED MASS TRANSPORT MODEL 71 5. COMPARISON AND DISCUSSION OF RESULTS 77 5.1 EFFECTS OF LATERAL DISPERSION 8 5.2 THE INITIAL DILUTION OF EFFLUENT 79 5.3 PEAK EFFLUENT CONCENTRATIONS 81 5.4 UPSTREAM EFFLUENT TRANSPORT 4 5.5 CHANNEL INTERACTIONS^ , 6 5.6 SUMMARY 86. SUMMARY AND CONCLUSIONS 92 REFERENCES. 98 -APPENDICES A. DERIVATION OF THE ONE-DIMENSIONAL MASS TRANSPORT EQUATIONS. . . .106 A.l GENERAL 10A.2 LONGITUDINAL ADVECTIVE TRANSPORT 107 A.3 LONGITUDINAL DISPERSIVE TRANSPORT 9 A. 4 SOURCE-SINK EFFECTS 113 A.5 THE ONE-DIMENSIONAL MASS TRANSPORT EQUATION 114 A. 6 THE LAGRANGIAN FORM OF THE ONE-DIMENSIONAL ' MASS TRANSPORT EQUATION 115 B. A DESCRIPTION OF THE FRASER RIVER ESTUARY 117 B. l GENERAL 11B.2 THE LOWER FRASER RIVER SYSTEM. . . 117 B.2.1 Channels of the River System .11B.2.2 Tides in the Strait of Georgia 129 B.2.3 River Discharges 12B.2.4 Saltwater Intrusion 134 vi APPENDIX Page C. NUMERICAL DISPERSION AND STABILITY 138 Cl SURFACE GEOMETRY OF PARTIAL DIFFERENTIAL EQUATIONS ..... 138 C.2 NUMERICAL DISPERSION 146 C.3 STABILITY . 154 C. 4 SUMMARY 159 D. DETAILS OF THE SOLUTION SCHEMES OF THE HYDRODYNAMIC AND MASS TRANSPORT EQUATIONS 161 D. l NUMERICAL SOLUTION OF THE HYDRODYNAMIC EQUATIONS 161 D.2 NUMERICAL SOLUTION OF THE TIDALLY VARYING MASS TRANSPORT EQUATION 164 D. 3 NUMERICAL SOLUTION OF THE TIDALLY AVERAGED 'MASS TRANSPORT EQUATION 7 E. ESTIMATION OF LATERAL DISPERSION 174 E. l EXISTING ESTIMATES OF LATERAL MIXING 17E. 2 VORTICITY ESTIMATE OF SECONDARY CURRENTS 180 P. ESTIMATION OF LONGITUDINAL DISPERSION 191 F. l GENERAL 19F.2 TIDALLY AVERAGED COEFFICIENTS OF LONGITUDINAL DISPERSION . . 192 F.3 TIME DEPENDENT LONGITUDINAL DISPERSION COEFFICIENTS 199 F.4 SENSITIVITY OF PREDICTED CONCENTRATIONS 201 F.5 SUMMARY 204 LIST OF TABLES Table Page 4.1 TIDALLY AVERAGED MASS BALANCE OF PREDICTED DISCHARGES FOR JANUARY 24, 1952 64 5.1 RATIO OF TIDALLY VARYING CONCENTRATION TO TIDALLY AVERAGED VALUE AT POINT OF DISCHARGE 85 B.l RIVER FLOW VOLUMES AND TIDAL PRISMS ON JANUARY 15 and JUNE 16, 1964 135 >E.l ESTIMATION OF COEFFICIENTS OF LATERAL DISPERSION 179 F.l ESTIMATED COEFFICIENTS OF LONGITUDINAL DISPERSION DUE TO LATERAL VELOCITY GRADIENTS 198 vii LIST OF FIGURES Figure Page 1.1 Effect of Variable Initial Dilution on Concentration 10 1.2 Effect of Multiple Dosing ;(Due to Flow Reversal) on Concentration 11 1.3 Movement of Effluent Through A Diverging Junction 14 3.1 The Hydrodynamic Estuary 36 3.2 Simplified Model Estuary of the Mass Transport Equations. . . 43 4.1 Equivalent Stage Response for Two Different Types of Cross-Section 50 4.2 Observed and Predicted Stages for June 16, 1964 52 4.3 Observed and Predicted Stages for June 16, 1964 53 4.4 Tidally Varying Stage and Discharge — Main Arm 55 4.5 Tidally Varying Stage and Discharge — Main Arm 56 4.6 Tidally Varying Stage and Discharge — Main Arm 57 4.7 Tidally Varying Stage and Discharge— Main Arm 58 4.8 Tidally Varying Stage and Discharge — North Arm 59 4.9 Tidally Varying Stage and Discharge — Middle Arm, Canoe Pass. 60 4.10 Tidally Varying Stage and Discharge — Pitt River 61 4.11 Tidally Varying Stage — Pitt Lake 62 4.12 Advection Of A Slug Load Down the Main Arm - Main Stem .... 66 4.13 Slug Inputs for Analytic and Predicted Dispersion Solutions. . 69 4.14 Dispersion of A Slug Load in the Main Stem 70 4.15 Influence of the Tidally Averaged Dispersion Coefficient on Predicted Tidally Averaged Concentrations .73 viii ix Figure Page 4.16 Maximum Upstream Excursion During Flow Reversal In Main Arm - Main Stem 75 4.17 Maximum Upstream Excursion During Flow Reversal In the North Arm and Pitt River 76 5.1 Initial Dilution at Station Nos. 10, 22, 50 and 102 80 5.2 Dispersion of A Concentration Spike 83 5.3 Predicted Concentrations in One Channel Caused by Effluent Discharge in Another Channel 87 5.4 Predicted Concentrations Due to Two Effluent Discharges in the Main Stem 88 A.l Elemental Cross-Sectional Slice of A River or Estuary. . . . 108 A. 2 Dispersive Effects of Vertical Velocity Gradients. ....... . 110 B. l The Drainage Basin of the Fraser River 118 B.2 The Fraser River From Hope to Vancouver 119 B.3 The Fraser River Delta 121 B.4 Network of Stations Used in the Numerical Solution of the Hydrodynamic Model . 124 B.5 Typical Channel Cross-Sections 125 B.6 Cross-Sectional Parameters of Main Arm - Main Stem 126 B.7 Cross-Sectional Parameters of the North Arm, Middle Arm and Canoe Pass 127 B.8 Cross-Sectional Parameters of. Pitt River and Pitt Lake. . . 128 B.9 Typical Tides at Steveston . 130 B.10 Local Low and High Tide:Envelopes 131 B.ll Tide Gauging Stations in the Fraser River Estuary 132 B.12 Mean Monthly Flows at Hope (1912-1970 Inclusive) 133 B.13 Salinity Profiles in the Main Arm on February 13-14, 1962 . 136 X Figure Page Cl Concentration Surfaces for the Advection of A Slug Load. . . 140 C.2 Time Rate of Concentration Along A Curve in the (x,t) Plane) 142 C.3 Dispersion of A Slug Load 144 C.4 Forward and Backward Time Differences 148 C.5 Upstream, Central and Downstream Space Differences 150 C. 6 Stable and Unstable Explicit Advective Schemes 158 D. l Explicit Finite Difference Grid of the Hydrodynamic Equations 163 D.2 Segments of Thomann's Tidally Averaged Model 168 D. 3 The Matrix A of Equation (D.8) for the Fraser River Estuary. 173 E. l Bends Along the Main Arm--: Main Stem 177 E.2 Tidally Varying Velocities .v 178 E.3 Assumed Linear Distribution of Secondary Velocities 181 E.4 Interaction of u and E, in A Straight River 184 z E. 5 Interaction of u and 5 Around A Bend 186 z F. l Lateral Velocity Profiles at Stations Nos. 14 and 15, Main Arm 194 F.2 Sensitivity of Tidally Varying Concentrations to the Coefficient of Longitudinal Dispersion 202 LIST OF SYMBOLS a : Height of estuary bed above level datum„ A : Cross-sectional area, A : Matrix of the effects of advection, dispersion and decay on the tidally averaged concentrations (Thomann's Solution). A^ ^+^: Average cross-sectional area between segments i and i+1 (Thomann's ' Solution) Aj : Cross-sectional area at point jAx at time nAt (hydrodynamic eguations) c i- Cross-sectionally averaged concentration b : Width of estuary C : Chezy's friction coefficient c : Concentration in segment i (Thomann's Solution) C : Matrix of segment concentrations (Thomann's Solution) Cj : Cross-sectionally averaged concentration at point jAx at time nAt. e : Coefficient of turbulent diffusion y E : Coefficient of longitudinal dispersion E^ : Coefficient of longitudinal dispersion due to the lateral velocity gradients of the steady velocity component Efc : Coefficient of longitudinal dispersion due to the lateral velocity gradients of the oscillatory velocity component E : Effective coefficient of longitudinal dispersion over a tidal c cycle due to lateral velocity gradients E : Coefficient of longitudinal dispersion due to the vertical velocity y gradients E^ : Coefficient of pseudo dispersion (numerical dispersion) Ei i+l: "^^HY averaged coefficient of longitudinal dispersion between ' segment i and i+1 (Thomann's Solution) xi Xll Gravitational acceleration Height of water surface above level datum Height of water surface above level datum at point jAx at time nAt Von Kantian's constant Rate of decay of substance in segment i (Thomann's Solution) Average length of segments i and i+l (Thomann's Solution) Manning's "n" (foot-second units) Freshwater discharge Radius of curvature of a bend Rate of production of substance due?to the i'th source-sink process Time Cross-sectionally averaged longitudinal velocity Magnitude of the steady component of the longitudinal velocity Amplitude of oscillatory component of the longitudinal velocity Effective longitudinal velocity over a tidal cycle Cross-sectionally averaged longitudinal velocity at point jAx at time nAt Shear velocity Shear velocity due to Shear velocity due to U Effective shear velocity over a tidal cycle Volume of segment i (Thomann's Solution) Cross-channel surface velocity due to secondary flows Rate of waste discharge into segment i (Thomann's Solution) Matrix of segment waste discharges (Thomann's Solution) Longitudinal distance along the estuary xiii y :. Vertical distance in the estuary z : Lateral distance in the estuary U*Y a : Factor given by a = ln{ —— } 11. C v a : Factor relating dispersion coefficients in depth and shear velocity, as in E = a y U# Ti^al exchange coefficient between segments i and i+1 (Thomann's Solution) Y : Standardized concentrations"'' predicted by the tidally averaged model Y^_v : Standardized concentrations^" predicted by the tidally varying model At : Time increment used in finite difference solutions Ax : Space increment used in finite difference solutions e : Coefficient of lateral dispersion : Coefficient of lateral dispersion due to the secondary flows associated with Uf e : Coefficient of lateral dispersion due to the secondary flows associated with U £c : Effective coefficient of lateral dispersion over a tidal cycle n : Relative depth, given by n = y/y £ : Vorticity generated in the lateral direction due to the vertical velocity gradients ?x : Longitudinal component of the lateral vorticity T : Circulation around a cross-section due to the longitudinal vorticity x "''The predicted concentrations are standardized by dividing the con centration determined from the mass of effluent and the freshwater flow per tidal cycle. ACKNOWLEDGEMENTS The author is very grateful to Dr. M. C. Quick for his guid ance, encouragement and criticisms during the development of this thesis. Mr. Donald 0. Hodgins, who developed and programmed the hydrodynamic model, is also thanked. Special thanks to go Mr. Richard Brun for draft ing the multitude of diagrams of this thesis, and to Miss Susan Aizenman for many pleasant typing and proofreading soirees. The latter part of this study was financially supported by the Westwater Research Centre at the University of British Columbia. The author expresses his appreciation for their assistance, with special thanks to Professor Irving K. Pox. Mr. F. A. Koch and Mr. G. S. Sheehan of the Westwater Research Centre are also thanked for their assistance during this study. Finally, the author gratefully acknowledges the patience and forbearance of his wife, Gloria, and son, Nicholas, during the course of this study. xiv INTRODUCTION Throughout history, centers of urban, agricultural and industrial development have commonly been located along rivers and estuaries. The principal uses of these surface water resources in early times were water supply and navigation; the river provided an accessible source of water for domestic and agricultural needs and a relatively cheap and easy means of bulk transport. Modern uses of surface water resources include naviga tion and water supply and waste disposal for domestic, agricultural and in dustrial purposes. In addition, rivers and estuaries provide habitat for wildlife, breeding and rearing areas for fish and shellfish and areas for general recreation. Associated with each use of surface water resources is a set of quantity and quality constraints that determine whether the water is satis factory for that particular use. The quality constraints consist., of the maximum allowable levels of various deleterious substances that may be pre sent in the water. Generally, the constraints for each use are different, and conflicts may arise when water is to.be used for multiple purposes. A common example is the conflict between the competing uses of waste disposal, wildlife habitat and recreation. It is recognized that the quantity and quality conflicts are related; ^ however, this thesis only considers quality aspects, and in particular, only those quality aspects that are determined One obvious way of improving water quality is by low flow augmen tation, as has been investigated by Worley et al.- [1965]. 1 2 by the concentration of a dissolved deleterious substance in the water. To investigate a situation of existing or potential water quality conflicts it is usual to'develop a "water-quality" model, or as it shall be referred to in this thesis, a mass-transport model. Such a model enables the concentration of deleterious substance to be predicted throughout the river or estuary, and is instrumental in assessing the effectiveness of possible control measures to improve water quality. The simplest type of estuarine mass-transport model is one-dimen sional (1-D) and only admits a variation of parameters and variables in the longitudinal direction. In a "tidally averaged" model, parameters and var iables are assigned their average values over a tidal cycle, whereas in a "tidally varying" model, they are allowed to vary throughout the tidal cycle. The temporal resolution of the tidally varying model is much finer than its tidally averaged counterpart, but considerably more effort is re quired in its development than for the latter. Thus, it seems relevant to enquire as to (1) whether the differences between the results of both models are significant; and (2) whether the extra effort involved in deve loping and applying the tidally varying model is justified by its finer resolution. In this thesis the significance of these differences is investi gated by applying both a tidally averaged and a tidally varying mass transport model to the Fraser River Estuary, a tidal estuary in the Province of British Columbia, Canada. The tidally averaged model was developed by Thomann [1963] and the tidally varying model was developed from first principles. Both mass transport models were developed as part of a larger interdisciplinary study 3 by the Westwater Research Center of the University of British Columbia to investigate the effects of possible patterns of future development on the water quality of the Fraser River Estuary. This thesis consists of six chapters. The tidally varying and tidally averaged forms of the one-dimensional mass transport equa^-tion are discussed in Chapter 1. The expected differencesLbetween the results of both models and the applicability of a one-dimensional model to the Fraser River Estuary are also considered,.there. In Chapter 2 the literature is reviewed to investigate the various ways of solving the one-dimensional mass transport equation. In order to solve the tidally varying form of the equation, it was necessary to develop a hydrodynamic model to predict the tidal variations in the longitudinal velocity and cross-sectional area along the estuary. The hydrodynamic model and the tidally varying and tidally averaged mass transport models applied to the Fraser River Estuary are described in Chapter 3. The verification of all three models is described in Chapter 4 and the results of applying both mass transport models to the Estuary are described and discussed in Chap ter 5. Finally, conclusions about the differences and applicability of both models are given in Chapter 6, There are six appendices to this thesis. The one-dimensional mass transport equation for unsteady non-uniform flow is derived in Appendix A. The advective and dispersive transport processes are discussed therein, as are'the assumptions in the derivation of the equation. Appendix B con sists of a description of the Fraser River Estuary. The various channels of the estuary, the freshwater flows, the tides and the salinity intrusion are 4 all described. When finite difference techniques are used to solve the mass transport equation,,there are problems with numerical dispersion and stability. In Appendix C both of these problems are seen to occur when a fixed grid is used to solve the equation, rather than a grid along the more fundamental characteristics of information propagation. Details of the finite difference schemes used to solve the hydrodynamic equations and both forms of the mass transport equation are given in Appendix D. In Appendix E, it is seen that existing estimates of lateral dispersion apparent ly underestimate the lateral mixing in the Fraser River Estuary. Secondary velocities are explained in terms of vorticity, and on the basis of very limited field data, the predicted secondary flows agree well with those ob served in the estuary. Finally, in Appendix F, estimates are made of the coefficients of longitudinal dispersion. An approximate method is given to allow for the variable contributions of the effects of vertical and lateral velocity gradients during the initial period before cross-sectional mixing is complete. CHAPTER 1 PRELIMINARY CONSIDERATIONS The one-dimensional mass transport equation is stated and its tidally varying and tidally averaged forms are briefly discussed. The problem of determining the parameters of both forms of the equation is considered, and expected differences between the results of a tidally varying and a tidally averaged mass transport model are described. The ability of a one-dimensional mass transport model to describe the mass transport processes in the Fraser River Estuary is also considered. 1.1 THE ONE-DIMENSIONAL MASS TRANSPORT EQUATION The one-dimensional equation for mass transport in unsteady non uniform flow in an estuary is obtained by taking a mass balance over an elemental cross-sectional slice of the estuary. Mass is transported through the slice by the mass transport processes of advection and dispersion, and these processes, together with any source-sink reactions that the substance undergoes, determine the concentration of the substance within the slice. As it is a one-dimensional equation, the dependent variable and the parameters are assigned their average cross-sectional values. The equation is derived in Appendix A and the assumptions in its derivation are discussed there. The one-? dimensional mass transport equation is given by |£ = - u|i +• ^{AEfe + ? S. (1.1 3t 9x A dx 9x . , 1 1=1 5 6 where c is the mean cross-sectional concentration of dissolved substance; u is the mean cross-sectional value of the longitudinal velocity; A is the cross-sectional area; E is the mean cross-sectional dispersion coefficient in the longitudinal direction; S,. is the rate of production per unit volume of water due to the itn source-sink process, it being assumed that there are n source-sink processes; x is the longitudinal distance; and t is the time. For the sake of brevity, the "mean" value of a parameter or variable is now taken to refer to its "mean cross-sectional" value. The high frequency turbu lent fluctuations are assumed to have been averaged out of u and c. The mean velocity u, the mean longitudinal dispersion coefficient E and the cross-sectional area A will be referred to as the parameters of Equation (1.1). These quantities are parameters in the sense that they are defined outside the equation by the cross-sectional geometry and hydraulics of the particular estuary. These three parameters can vary with both x and t Thus, Eguation (1.1) is seen to be a second-order linear partial differential equation with variable coefficients. The three terms on the right-hand side of the equation will be referred to as the advective, the dispersive and the source-sink terms respectively. 7 Equation (1.1) was derived for the general case of unsteady non uniform flow and is the basis of all one-dimensional mass transport models of rivers and estuaries. In an estuary, the rise and fall of the tide causes temporal variations in the parameters u, A and E. In a tidally varying model, this temporal variation in the parameters is taken into account as Equation (1.1) is solved to determine the mean concentration along the estuary during the tidal cycle. In a tidally averaged model, the parameters are assigned their average values over a tidal cycle, and the mean concentration along the estuary is determined over periods of a tidal cycle. Averaging Equation (1.1) over a tidal cycle does not alter the form of the equation, but merely changes the interpretation of the^dependent variable and the parameters. For example, u becomes the tidally averaged velocity and is determined by the freshwater discharge through the tidally averaged area A. it should be noted that the concentration predicted by a tidally averaged model is not the "tidally averaged" concentration — this can only be determined by aver aging the results of a tidally varying model over the tidal cycle. The majority of tidally averaged models supposedly predict the mean concentration along the estuary at times of slack-water. The reason for this is the ease of sampling at times of slackwater, and is discussed further in Section 2.-1-: 1.2 DETERMINATION OF PARAMETERS The flow field of an estuary consists of an unsteady oscillatory component due to the tide superimposed on a steady component due to freshwater inflow. Before a solution can be obtained to the tidally varying form of Equation (1.1), it is necessary to determine the tidally induced temporal 8 variation in the parameters u, A and E. A one-dimensional hydrodynamic model was developed to determine the temporal variation in u and A. In this model, the equations of motion and continuity were applied to the water mass of the estuary and solved throughout the tidal cycle. The hydro-dynamic model is described in detail in Chapter 3. The temporal variation in E during the tidal cycle was related to the temporal variation in u, as is discussed in Appendix F. In a tidally averaged model, the parameters u, A and E are aver aged over a tidal cycle. This reduces the unsteady tidal flow field to a steady freshwater flow field, and consequently u can be determined from the equation of continuity alone (the equation of continuity is applied to the freshwater discharge through the tidally averaged area A). The tidally averaged dispersion coefficient includes the effects of upstream advection on the flood tide, as"is discussed in Chapter 4. With a tidally varying model, mean concentrations are determined throughout the tidal cycle, whereas with a tidally averaged model, mean concentrations are determined over periods of a tidal cycle. In effect, the temporal resolution of the tidally varying model is much higher than that of a tidally averaged model. Because of their unsteady nature, the estimation of the tidally varying parameters u and A involves significantly more effort than the estimation of their steady tidally averaged counter parts. This is apparent from the discussion of Chapters 3 and 4. Thus the greater temporal accuracy of the tidally varying model is offset by the greater effort required to estimate the tidally varying parameters. 9 1.3 THE EFFECTS OF THE TIDE ON MASS TRANSPORT PROCESSES In an estuary, the tidal rise and fall of the water surface causes temporal variations in the longitudinal flow. Because of this temp oral variation of flow, the initial dilution of a discharged effluent also varies during the tidal cycle. This is illustrated in Figure 1.1 for the steady discharge of effluent into an estuary during the flood tide. The effect of the flooding tide is to continuously reduce the seaward flow of water past the effluent outfall, and thus the slug of water dosed in any time increment becomes increasingly smaller. However, since the effluent discharge is steady, the same mass of effluent is added to each of the dosed slugs and this gives rise to the spatial distribution of cross-sectionally averaged concentration shown in Figure 1.1. As the tide continues to flood, flow reversal will occur at the effluent outfall and previously dosed slugs will move upstream past the outfall and be dosed again, as is illustrated in Figure 1.2. On the following ebb tide, upstream slugs will move down stream past the outfall and be dosed yet again, as is also shown in Figure 1.2. A tidally varying mass transport model can account for the effects of variation in initial dilution and multiple dosing, whereas a tidally aver aged model cannot. Essentially, the effect of the tidal variation in flow is to in troduce "spikes" into the concentration profile along the estuary (as is apparent from Figures 1.1 and 1.2). The concentration of these spikes is then reduced by the dispersive transport process. Thus, differences between the results from a tidally varying and tidally averaged mass transport model will depend on the magnitude of the tidally varying dispersion process. Another difference between the results of both models is that a tidally vary ing model correctly accounts for the effects of upstream advection on the 10 Steady effluent discharge t = 0 Rising Tide ESTUARY SEA t = At Cl f t = 2At C2 Cl t = 3 At C3 C2 Cl 1 C4 t = 4At C4 C3 C2 Cl u Time varying dilution C3^Time average dilution Figure 1.1 Effect of Variable Initial Dilution on Concentration 11 Steady effluent discharge t = 4At i Rising Tide C4 C3 C2 Cl SEA u = 0 f Rising Tide t = 5At C4 C5 C3 C2 Cl Falling Tide t = 6At C|C3 C2 Cl 1 WC4+C5 + C6 Effect of multiple dosing Time average concentration m////////////////A ci Figure 1.2 Effect of Multiple Dosing (Due to Flow Reversal) On Concentration 12 flood tide, whereas a tidally averaged model can only simulate this upstream transport via the tidally averaged dispersion process. 1.4 ACCURACY OF A ONE-DIMENSIONAL MODEL Having briefly discussed the tidally varying and tidally averaged forms of the mass transport equation, the problem of how well a one-dimension al model describes the mass transport processes in the Fraser River Estuary is now considered. When effluent is discharged into a river or estuary, the time required for complete cross-sectional mixing to occur depends on the coefficient of lateral dispersion, the width of the river or estuary and the position of the effluent outfall [Ward, 1973]. In the initial period before cross-sectional mixing is complete there are significant lateral concentra tion gradients across the river or estuary. A one-dimensional mass trans port model predicts the cross-sectionally averaged concentrations and does not "see" these lateral gradients. Thus, in the initial period before cross-sectional mixing is complete, the predicted concentrations will underestimate the peak lateral values in the river or estuary. Also, during this initial period the dispersion of effluent results in a skewed distribution of concen tration along the river or estuary rather than the Gaussian distribution pre dicted by the tidally varying mass transport model (see Appendix A). When the cross-sectional mixing is complete, longitudinal gradients dominate the transport processes, and a one-dimensional model will provide a good descrip tion of the actual concentration profile along the river or estuary. 13 The Fraser River estuary consists of the seven principal channels described in Appendix B. The junctions of these various channels influence how effluent is advected through the estuary. If cross-sectional mixing is not complete at a diverging junction, the bulk of the effluent load may be advected down one of the channels as shown in Figure 1.3. A one-dimensional mass transport model cannot reproduce this behaviour at the junction. The situation shown in Figure 1.3 is further complicated by the effects of tidal flow reversal and the presence of secondary flows at the junction. On the flood tide some effluent is carried back upstream past the junction, and additional cross-sectional mixing occurs before the effluent is advected back through the junction on the ebb tide. If secondary flows are present, the situation may arise where effluent is released from one bank upstream from a diverging junction and is then advected down the opposite channel, as is shown in Figure 1.3. (On the ebb tide, marked secondary flows of this nature are commonly observed at the Main Arm-North Arm junction shown in Figure B.3). Thus, in the period before cross-sectional mixing is complete, the movement of effluent through the various junctions of the Fraser River Estuary is a complex two-dimensional process that also varies during the tidal cycle. To reproduce these effects, a one-dimensional model would have to be modified at the junctions and such modification may require con siderable field data. It is noted that the effects of tidal flow reversal and secondary flow will enhance the cross-sectional mixing at junctions. The ability of a one-dimensional mass transport model to describe the concentration profile along the various channels of the Fraser River Estuary essentially depends on the time required for cross-sectional mixing to occur. In the initial period before cross-sectional mixing is complete, Movement of Figure.1.3 Effluent Through A Diverging Junction 15 the dispersion process is skewed rather than Gaussian, the peak lateral concentration is significantly greater than the predicted cross-sectionally averaged value and the movement of effluent through the junctions may not be according to the simple mass balance of the one-dimensional equation, the latter two effects probably being the most important. The time required for 80 per cent cross-sectional mixing to occur has been estimated in Appen dix E. The techniques of Fischer [1969a] and Ward [1972] gave an estimate of 55 hours, whereas calculations based on vorticity considerations gave an estimate of five hours, the bulk of the cross-sectional mixing being due to the influence of secondary flows. It is noted that this last value has not been confirmed by field experiments. Because of tlje highly asymmetrical nature of the tides of the Fraser River Estuary, a more realistic estimate of the time of cross-sectional mixing in the lower reaches of the estuary is probably 1-2 tidal cycles. All of these aspects of lateral mixing are discussed in Appendix E. Unfortunately, time and expense have precluded using dye studies to measure the actual rate of cross-sectional mixing in the estuary. To sum up, in this study it is recognized that lateral mixing and the complex two and three-dimensional flow characteristics through the junctions can have a considerable influence on the concentrations along the various channels of the Estuary. These aspects are considered', and their influence is partially assessed, but the main thrust of the study has been to obtain accurate solutions to the tidally varying and tidally averaged forms of the one-dimensional mass transport equation so that the influence of the tides on the predicted concentrations can be assessed. To some extent, 16 Che two- and three-dimensional effects can be estimated as further modifi cations of the one-dimensional results, as is discussed further in Chapter 5. CHAPTER 2 LITERATURE REVIEW The literature is reviewed to investigate the various solutions and applications of one-dimensional mass transport models. Solutions to the mass transport equation are classified into the categories of analy tical solutions, numerical solutions, physical and analogue model solutions and stochastic solutions. The advantages and disadvantages of each class of solution are also discussed. 2.1 ANALYTICAL SOLUTIONS The one-dimensional mass transport equation for unsteady non uniform flow in a river or estuary is given by |£ = -up- + I^-fAE^} + E S., (2.1 3t 3x A 9x dX . . l' 1=1 where the terms are as defined in Section 1.2. Equation (2.1) is a second order linear partial differential equation with variable coefficients (the source-sink terms are generally linear,- and u, A and E are dependent on x and t). As such, no completely general analytical solution exists, but solutions have been obtained under a number of simplifying assumptions. Before reviewing various solutions, the so-called "slackwater" concentra tions predicted by the tidally averaged models are discussed. The majority of the analytical solutions discussed here are for the steady-state response (dc/dt = 0) of Equation (2.1). This simplification 17 18 reduces the partial differential equation to an ordinary differential equation in x alone. In determining this steady state".response, the freshwater discharge and effluent inputs are assumed to remain steady or constant for a period of time equal to the residence time of the estuary. Further, any steady-state solution of Equation (2.1) is a tidally aver aged solution. (The equation is now independent of time and cannot re solve the within tide temporal fluctuations in c, u, A and E). However, the steady state solution to Equation (2.1) does not represent the tidally averaged concentration profile along the real estuary. In fact, it repre sents the concentration profile along a model estuary that is tideless and has a high degree of longitudinal dispersion. (The tidally averaged E is generally much higher than the tidally varying E, as is discussed in Sec tion 3.3). The tidally averaged concentration profile can only be obtained from the results of a tidally varying mass transport model. In the majority of tidally averaged mass transport models, the steady state solution to Equation (2.1) is assumed to represent the concen tration profile along the estuary at times of slackwater. In other words, the tidally averaged mass transport processes are used to determine the concentration profile along the estuary at a particular phase of the tide. The reason for working with slackwater profiles is the relative ease of sampling tracers in the estuary at these times [;0'Connor, 1960] . Any corres pondence between the steady-state solution of Equation (2.1) and the slack-water concentration profile along the real estuary is because the tidal effects are small, and under these conditions the steady-state solution is an adequate representation of the concentration profile in the real estuary 19 at any phase of the tide, or because the steady state solution is "forced" to conform to the slackwater profile by "adjusting" the parameters of Eguation (2.1). It is interesting to note that Preddy and Webber [1963] developed a tidally averaged model that predicts tidally averaged concentra tions rather than slackwater values. This model is described in Section 2.2.2. In a tidally varying situation, quasi-steady-state conditions are said to be achieved when the variable of interest (for example, c) undergoes a cyclical repetition of the same values from tide cycle to tide cycle. Thus, in a tidally varying model one speaks of a quasi-steady-state response, rather than a steady-state response as with a tidally averaged model. Essentially, the analytical solutions of Equation (2.1) entail simplifications in which various terms of the equation are ignored (for example, the dispersion term) and the parameters of the u, A and E are repre sented as simple functions of x and possibly t. Effluent inputs are treated as boundary conditions, and solutions have been obtained for both point and distributed effluent inputs [O'Connor, 1965] . Steady-state solutions to Equation (2.1) have been obtained for steady effluent discharge into (1) steady uniform flows [O'Connor, 1960, 1962]; (2) steady non-uniform flows where the area is a simple function of x [O'Connor, 1965, 1967]; and (3) steady non-uniform flows where the freshwater discharge varies exponentially with x due to land run-off [O'Connor, 1967]. Transient solutions have been obtained for steady effluent discharge into (1) non-steady freshwater flows during the recession limb of the hydrograph (which was approxi-20 mated as negative exponential in time) [O'Connor, 1967]; and (2) steady flows where the diurnal photosynthetic production of DO is taken into account (it was approximated as a half sine wave) [O'Connor and Di Toro, 1970]. To obtain these transient solutions, it was necessary to ignore dispersion. Kent [1960] used the method of separation of variables to ob tain the general solution to Equation (2.1) for a slug input into steady uniform flows, and the particular solutions for a slug input into steady non-uniform flows (u, A and E were assumed to be linear in x). Di Toro and O'Connor [1968] obtained the transient solution for steady effluent input into unsteady non-uniform flows where the variation in cross-sectional area could be separated into independent functions of x and t (dispersion was ignored). Li [1962] used the method of characteristics to obtain a solution for non-steady effluent discharge (sinusoidal variation) into steady uni form flow-, (dispersion was ignored) . He later used the method of perturba tions to obtain a solution for the same case with dispersion included [Li, 1972] . Holley [T969b] transformed the mass transport equation for BOD into the elementary diffusion equation. (By travelling with the water mass, the only transport process an observer "sees" is dispersion — see Section A.6). He took the standard solution for a slug input and used the convolution method to obtain the solution for the continuous (but not neces sarily steady) effluent discharge into unsteady uniform flow. Bennet [1971] used the convolution method to obtain the solution for the complete BOD-DO 21 system for the same effluent discharge and flow conditions as Holley. 2.2 NUMERICAL SOLUTIONS 2.2.1 Finite Difference Solutions. Finite difference methods for the solution of partial differential eguations can be classified into (1) fixed mesh methods in which the solution is obtained at fixed predetermined points in a rectangular mesh of time and distance; (2) characteristic methods, in which the solution is obtained at mesh points along the characteristic curve(s) in the time-distance planets), the position of the mesh points being determined as the solution progresses; and (3) combined methods, in which the solution is followed along the characteristic curve(s) and then extrapolated back onto a fixed mesh of uniformly spaced points [Amein, 1966]. In characteristic methods, the mesh points are generally non-uniformly spaced in time or distance, and consequently the book-keeping of results is somewhat untidy. Combined methods simplify this bookkeeping by extrapolating the results back onto a fixed uniformly spaced mesh. Finite differenceemethods are discussed in some detail in Section C.2. The problems of stability and convergence are considered, and it is seen that implicit finite difference schemes are generally unconditionally stable, whereas explicit schemes are at most conditionally stable. The stability requirements of explicit schemes are discussed in detail in Section C.3, and are seen to impose limits on the relative size of Ax and At, the grid spacing. Generally, fixed mesh finite difference schemes do not simulate the advective transport process correctly, and result in an additional dis persive process being superimposed on the actual advective and dispersive 22 processes occurring in the rivertor estuary. This so-called numerical dispersion is discussed in detail in Section: c.2 and is seen to be elimi nated by using characteristic finite difference methods. Harleman et al. [1968] developed a one-dimensional, tidally varying mass transport model of the tidal portion of the Potomac River. A fixed mesh, implicit scheme with central differences was used to solve the mass transport equation. (The various types of space and time differ ences are described in Section C.2). Tidal velocities were calculated from discharge and tidal records, and the Taylor dispersion equation [Taylor, 1954], modified for the effects of vertical velocity gradients, was used to determine the longitudinal dispersion coefficients. The finite differ ences?' solution simulated the results of dye studies in the estuary reason ably well. However, Prych and Chidley [1969] showed that the numerical dispersion in their finite difference scheme was approximately 30 times greater than the modified Taylor dispersion, and was of the order of the same magnitude as the actual dispersion occurring in the estuary. (The modified Taylor equation neglects the effects of transverse velocity grad ients and grossly underestimates the dispersion coefficient for streams and estuaries, as Fischer [1966b] has demonstrated). The agreement between the dye results and the finite difference solution was apparently only fortui tous. This example illustrates the significance of numerical dispersion. Bella and Dobbins [1968] investigated fixed mesh finite differ ence solutions to the one-dimensional mass transport equation describing the BOD - DO system. They applied a tidally varying finite difference model to a hypothetical estuary and compared the results with an analytical tidally averaged solution. 23 Dornhelm and Woolhiser [1968] obtained a solution to the one-dimensional, tidally varying mass transport model for conservative sub stances. They used a fixed mesh implicit scheme with central differences. A hydrodynamic model, solved by the same finite difference scheme, was used to determine the tidally varying parameters. The hydrodynamic model was verified against a steady-state analytical solution, but it exhibited instabilities when applied to the Delaware estuary. (According to a linear stability analysis, their implicit scheme was unconditionally stable. How ever, the non-linear nature of the hydrodynamic equations may require more stringent stability conditions, as is discussed in SectionD.l). The mass transport model was verified for steady uniform flow and was applied in tidally varying form to a hypothetical estuary. To overcome the problems of numerical dispersion, Gardiner et al. [1964] used the method of characteristics to solve the two-dimensional mass transport equation describing the movement of a solvent through sand saturated with oil. A combined finite difference method was used, with a mesh of uniformly distributed moving points superimposed on a fixed two-dimensional space mesh. At the beginning of a time step, the moving points are assigned the concentration of the nearest fixed grid point, and then ad vected along the characteristics according to the equation analagous to Equation (A.9). After advection, the moving points are assigned to the closest grid point to determine the concentration there. An explicit form of the equation analagous to Equation (A.8) is then used to disperse the concentration over the fixed grid points. The complete process is then repeated for the next time step. Pinder and Cooper [1970] used the same 24 technique to calculate the transient positions of the salt water front in a coastal aguifier. Di Toro [1969] recognized that-the method of characteristics re duced the mass transport equation for the BOD-DO system into two ordinary differential equations. (He ignored dispersion, and so the equation anal-agous to Equation (A.8) was a "true" ordinary differential equation). He noted that the numerical solution of ordinary differential equations is more exact than for partial differential equations, and that the accuracy can be controlled by predictor-corrector methods. 2.2.2 "Box Model" Solutions. Another type of mass transport model is the so-called "box model" in which the estuary is divided into finite segments or boxes. Each segment is assumed completely mixed and the concentration of substance in a segment is determined by the discharge of substance into the segment, the advective and mixing processes occuring between adjacent segments, and any source-sink effects that the substance undergoes. Callaway [1971] noted the box models are conceptually very similar to finite difference models, as in effect, the latter segment the estuary into well mixed boxes centred around the grid points. Another similarity is that the box model representation of the estuary results in a set of difference equations relating the concentration in any box to the concentration in neighboring boxes (similar to implicit difference schemes). Thomann [1963], [1965], developed a tidally averaged box model which was applied to the Delaware estuary. The estuary was segmented into 30 boxes whose length varied from two to four miles. All hydraulic parameters 25 were assumed steady, but the concentration within each box was allowed to vary with time so that the effect of time varying inputs could be studied. The concentration in each box is determined from a set of simultaneous difference-differential equations. For' steady-state inputs, the temporal derivatives are set equal to zero and the system of equations reduces to a set of difference equations similar to those of an implicit finite differ ence scheme. Thomann's model is described in more detail in Sections 3.3 and D.3. Pence et ai. [1968] extended the model to account for time varying fresh water inflows(tidal parameters are still averaged over a tidal cycle). Like most tidally averaged models, Thomann's model is a "slack water" model with concentrations being determined at succeeding slack waters. Preddy and Webber [1963] developed a tidally averaged box model of the Thames estuary to predict the true tidally averaged value of the DO concentration. In their model, the boxes were two miles long and the tidal displacement was six miles to either side of a box. Thus, the tidally averaged concentration within a box depends on the concentrations within the three boxes upstream and downstream of it. This model is one of the few that determine tidally averaged concentrations rather than slack water values. 2.2.3 Finite Element Solutions. Finite element methods have been used to obtain solutions to the one-dimensional mass transport equation. Examples include Price et al. [1968], Guymon [1970], Pinder and Friend [1972] and Nalluswami et al. [1972]. In finite element solutions, the estuary is segmented and the concentration profile is approximated as a series of func tions. 26 In. some respects the method is similar to implicit finite differ ence methods and box models; in the latter two methods the estuary is also segmented, but the concentration profile is approximated as a point value rather than a function in each segment. Further, the solution of all three methods is given by a set of simultaneous equations governed by a square banded matrix of bandwidth three (see Section D.3 for Thomann's version of this matrix). Price et al. [1968] compared finite element, finite differ ence and method of characteristics solutions for the case of a constant discharge into steady uniform flow. The finite element method was found to be more accurate than the finite difference approximations, but this is only to be expected as this technique can follow variation within a segment, whereas finite difference solutions cannot. The finite element method was faster and more accurate than the method of characteristics of Gardiner et at. [1964] However, the accuracy of the latter technique was found to be adequate, and any inaccuracies probably arise from the procedure of extrapolating the moving points back onto the fixed grid points. Apparently, there are no problems with stability and numerical dispersion in finite element methods [Price et ai., 1968]. In closing this section, it is mentioned that two-dimensional finite difference finite element and box models have been developed. Fischer [1970], Oster et al. [1970] and Leendertse [1971a] are among those who have obtained finite difference solutions for two-dimensional mass transport in the horizontal plane, whereas the box model of Pritchard [1969] was developed to describe two-dimensional mass transport in the vertical plane. 27 2.3 PHYSICAL AND ANALOGUE MODEL SOLUTIONS 2.3.1 Physical Model Solutions. Physical models of estuaries were origi nally developed to investigate sediment erosion and deposition in tidal water ways. The use of these models to investigate mass transport processes in the prototype estuary was a natural development, and examples include O'Connor [1962], [1965]; Diachishin [1963]; O'Connell and Walter [1963] and Lager and Tchobanoglous [1968]. In these model investigations of the mass transport process, a quantity of dye is introduced at the point being investigated and its distribution throughout the model estuary is recorded over succeeding tidal cycles. Physical modelling involves the scaled-down reproduction of the more important processes that affect the parameter being modelled. For mass transport in estuaries, these processes are advection, dispersion and source-sink effects. Physical models can satisfactorily reproduce one-dimensional advection and the two-dimensional density dependent salinity intrusion pro cess [Simmons, 1960; Harleman, 1965]. However, the use of distorted cross- . sectional space scales in the model results in distorted cross-sectional velocity distributions. From the discussion of Section A.3, it is apparent that this will result in an incorrect reproduction of-the! dispersion process. This has been discussed by Harleman [1965], Harlemen-et al. [1968], and Fischer and Holley [1971]. Even in physical models with undistorted cross-sectional space scales, the dispersion process may not be reproduced correctly. In this case, the reproduction of the dispersion process depends on the ratio of the time of cross-sectional mixing to the tidal period in the prototype estuary [Fischer and Holley, 1971]. 28 No attempt is made to reproduce any source-sink reactions in model studies, a "conservative" dye being used for the tests. However, the dye does adsorb onto surfaces, and this must be allowed for in inter preting test results. O'Connell and Walter [1963] developed a method to account for a first-order decay reaction. 2.3.2 Analogue Model Solutions. Electrical analogue models have been used to obtain solutions to the mass transport equation [Rennerfelt, 1963 and Leeds, 1967] and the hydrodynamic equations [Harder, 1971]. The estuary is divided into segments, and the space derivatives of the mass transport equation are approximated as finite differences over those seg ments. This reduces Equation (2.1) to a set of simultaneous difference-differential equations (the time derivative is continuous) similar to those of Thomann [1963]. The electrical analogue of this system of equations can then be constructed in the form of a so-called ladder network [Leeds and Bybee, 1967]. Analogue models can determine the steady- state and transient solutions for either constant or sinusoidally varying effluent discharge into steady non-uniform flows [Leeds, 1967 and Leeds and Bybee, 1967], As the flow is assumed steady, these models are tidally averaged, but it may be possible to develop an analogue for sinusoidally varying flows. Depending on the finite difference approximations used for the spatial derivatives, these models may also suffer from numerical dispersion, as was recognized by Bella [1968]. For example, the analogue of Leeds and Bybee [1967] used central differences to approximate the space derivatives of Equation (2.1). From the discussion of Section C.2, it is apparent that this will result in numerical dispersion. 29 2.4 STOCHASTIC SOLUTIONS Before reviewing stochastic solutions of the mass transport equa tion, it is necessary to define several terms. A process is any phenomena that undergoes changes with respect to time, an example being daily river flow. If the chance of occurrence is taken into account, a process is said to be stochastic or probabalistic. A probabilistic process is time-independent and the variables are considered pure-random. A stochastic process is time-dependent and the variables may be pure-random or non-pure-random. If non-pure-random, the process is composed of a deterministic and a pure-random component. If the probability distribution of the random variable remains constant throughout the process, the process is said to be stationary, other wise it is non-stationary. Non-stationary stochastic processes are very complicated mathematically, and in order of increasing simplification they are treated as stationary, probabalistic, and deterministic [Chow, 1964]. Stochastic mass transport models attempt to predict both the mean concentration profile and the variation around the mean. This varia tion is due to the stochastic nature ofC.thV underlying processes that deter mine the concentration, namely the variations that occur in the freshwater and effluent discharges, effluent concentrations, etc. during the period of analysis. Deterministic models use the mean value of each of these stochastic variables to determine the mean concentration profile. Diaschishin [1963] assumed tidal mixing to be a pure-random pro cess that can be characterized by a mixing length. On this basis he used a random walk formulation to obtain a tidally averaged solution for waste dis posal in purely tidal waters. He allowed for the effect of mixing in the \ver-30 tical and lateral directions, as well as in the longitudinal direction. The solution of Equation (2.1) for a slug input into steady uniform flow is a Gaussian distribution. Harris [1962, 1963] assumed this distribution to result from a pure-random dispersion process, and on this basis he obtained maximum likelihood estimates of u and E. He used the convolution method to obtain the solution for a continuous re lease. However, in the initial non-Fickian period, the cross-sectionally averaged concentration profile is skewed and not Gaussian (see Section A.3), as Fischer [1966a] noted. Loucks and Lynn [1966] obtained the transient and steady state probability distributions of the DO concentrations at the point of minimum DO in a stream (dispersion was ignored). The sequence of daily streamflows was assumed to be a first-order Markov process (a stochastic non-pure-random process), and the sewage flows, ultimate BOD and source-sink parameters were assumed to be dependent on the daily streamflows (via conditional probabili ties) . Essentially, the technique is as follows. Given a streamflow, the set of sewage flows, source-sink parameters, etc. that result in the minimum DO falling below some prescribed value is determined (by a modified Streeter-Phelps equation). The probability of this set of events occurring is then determined, the process being repeated for each discrete streamflow variate, and the probabilities summed to give the total probability of the minimum DO being less than the prescribed value. Loucks [1967] used a similar technique to investigate the effect of various treatment plant operating policies on the distribution of minimum DO in a stream. Thayer and Krutchkoff [1967] obtained a stochastic solution for the concentration profiles of BOD and DO in a stream (dispersion was ignored, and 31 the streamflow was assumed uniform and steady). The basic processes of decay, reaeration, etc. occurring in the modified Streeter-Phelps eguation were treated as stochastic, but freshwater and effluent discharges were assumed to be deterministic. Concentration values were divided into dis crete units of size A, and the probability of a change of A occurring in the concentration during a small time interval was assumed to follow a Pois son distribution. This allowed them to lump all of the stochastic variation into the parameter A (the sum of a number of independent Poisson processes is itself Poisson), which was determined by fitting the predicted variance in the model estuary to the observed variance in the real estuary. Essentially, they used the same set of results to both estimate and verify the model. Under these conditions, any model will reproduce the observed results, irres pective of whether the underlying processes are correctly modelled. Whereas other stochastic solutions use the explicit variation in the underlying stochastic processes to predict the variation in the concentration profiles, Thayer and Krutchkoff based their solution on the implicit variation in the stochastic processes (the actual variation is not even measured). Their predicted mean concentration value is simply the:'solution to the (determinis tic) modified Streeter-Phelps equation. As such, their technique is no better than any deterministic model that recognizes and incorporates stochastic vari ation through measurements of field values. Custer and Krutchkoff [1969] used a random walk formulation to extend the technique to include dispersion (the flow was assumed uniform but unsteady), and Schofield and Krutchkoff [1972] further extended the technique to account for variable stochastic parameters along the estuary. 32 Koivo and Phillips [1971] regarded the measured BOD and DO concen tration profiles along a stream to be corrupted by "noise" due to the stochas tic nature of the underlying processes. They developed a method based on non linear regression analysis to determine the values of the stochastic para meters that gave the best overall fit to the observed concentration profiles. 2.5 SUMMARY Most analytical solutions to the one-dimensional mass transport equation are tidally averaged and may be useful for investigating effluent discharge into estuaries with simple geometries. (When the number of segments with different geometries is greater than four or five, the matching of bound ary conditions at their ends becomes cumbersome [O'Connor et al, , 1968]) . In tidally averaged solutions, upstream transport is by dispersion, and thus a tidally averaged model may not adequately reproduce the interaction of two or more effluent outfalls on the flood tide. Numerical solutions provide a means of solving the one-dimensional tidally varying mass transport equation, but the problems of stability and numerical dispersion must be considered. Of the various numerical methods, the method of characteristics, methods with higher order approximations and finite element methods are satisfactory with respect to accuracy of solution. The method of characteristics has the additional advantage of directly simula ting the advective process. Undistorted physical models may be useful in resolving complex three-dimensional aspects of the flow-field at important or sensitive sections of the estuary, but because of their inaccurate and inconsistant reproduction 33 of the dispersion process (in both distorted and undistorted models), at best they are only an adjunct to a mathematical form of solution. As deve loped, analogue models are tidally averaged and suffer from numerical dis persion, although it may be possible to overcome both these limitations. It is not clear whether the various stochastic solutions attempt to model "true" stochastic variation, or a combination of stochastic and cross-sectional variation. If the cross-sectional variation is greater than the actual stochastic variation, it may be more meaningful to model the lateral component of the dispersion process with a two-dimensional deterministic model. CHAPTER 3 A DESCRIPTION OF THE HYDRODYNAMIC AND MASS TRANSPORT MODELS 3.1 THE HYDRODYNAMIC MODEL To solve the one-dimensional tidally varying mass transport equation, it is necessary to know both the spatial and temporal variation in the parameters u, A and E of the eguation, as is discussed in Section 1.3. A one-dimensional hydrodynamic model was developed to predict the spatial and temporal variation in the parameters u and A. (The temporal and spatial variation of the parameter E is discussed in Appendix F. In the hydrodynamic model, the equations of motion and continuity were applied to the water mass of the estuary, and solved throughout the tidal cycle. River Estuary, the equations of motion and continuity are given by [Dronkers, 1969] 3.1.1 The Hydrodynamic Eguations. As applied to the Fraser 9u 3t dh u u (3.1) + u- •g. (3.2) where u is the mean longitudinal velocity; h is the height of the water surface above an arbitrary level datum; y is the mean cross-sectional water depth; 34 35 A is the cross-sectional area; b is the cross-sectional width; g is the local gravitational acceleration; and C is Chezy's friction factor. Several of these terms are illustrated in Figure 3.1. Note that in the hydrodynamic equations, x increases in the upstream direction, whereas in the mass transport equation, x increases in the downstream direction. The two terms on the left-hand side of Equation (3.1) are the local and convective (or Bernoulli) accelerations respectively. The two terms on the right-hand side of the equation are the forces causing these accelerations, the net pressure force due to the slope of the water surface and the friction force respectively. The first term of Equation (3.2) is the net outflow from an-;elemental cross-sectional slice of the estuary, while the second term represents the accompanying change in storage within the slice. The assumptions made in deriving equations (3.1) and (3.2) are listed in Dronkers [1964]. The most important assumption is that the storage width is equal to the advective width. The variation of advective and stor age widths along the major channels of the estuary is shown in Figures B.6 to B.8. The difference between the advective and storage widths is exagger ated since these values were determined from the cross-sectional areas below local low water. (In the lower reaches of the estuary, the tidally averaged depth is some 10 feet above local low water). Consequently, the assumption of equal advective and storage widths seems reasonable for the Fraser River Estuary. (This is also discussed in Section 4.1 with regard to the verifi-36 Figure 3.1 The Hydrodynamic Estuary 37 cation of the hydrodynamic model). The one-dimensional nature of the equations is justified since the various channels of the estuary are much longer than they are wide (see Appendix B). In fact, Callaway [1971] has classified this type of estuary as a "tidal river." In deriving Equa tions (3.1) and (3.2) the presence of the salt wedge has been ignored be cause of the complicated nature of its dynamics (see Appendix B). Accord ing to Odd [1971], the Chezy formula should adequately represent frictional effects in fast flowing, well mixed estuaries. The Fraser River estuary is fast flowing, but the Chezy formula may not be accurate in the highly stra tified lower reaches. From Figure 3.1, it is seen that h = a + y (3.3) where a is the height of the river bottom above the same level datum that h is referred to. Thus, Equations (3.1) and (3.2) are seen to be a pair of coupled partial differential equations, the dependent variables being u and either y or h, and the independent variables being x and t. The width b, Chezy's friction factor C, the area A and the factor a all appear in the equations as parameters. In deriving these equations, it has been assumed that the dependent variables and parameters vary in a continuous manner along the estuary. It is noted that both equations contain variable coeffi cients and that Equation (3.1) contains non-linear terms. The fixed mesh, explicit finite difference method of Dronkers [1969J was used to obtain a numerical solution to the hydrodynamic equations. De tails of the method are given in Appendix D and the fixed mesh of "stations" used in the finite difference solution is shown in Figure B.4. Because the 38 solution scheme of the hydrodynamic equations is explicit, stability re quirements determine the relative size of Ax and At. In solving the hydrodynamic eguations, Ax was set equal to 5,000 feet, except in the deeper waters of Pitt Lake where Ax was increased to 15,000 feet for stability reasons. It was found that with this space grid, a At of 90 seconds was satisfactory as regards stability. The stability requirements are discussed in detail in Section C.3. 3.1.2 Assumed Quasi-Steady Hydraulic Conditions. For effluent discharge into an estuary, the quasi-steady state response of an estuary is usually a condition of interest. .Tidally. varying mass transport models usually require many tidal cycles to achieve quasi-steady-state conditions, whereas hydrodynamic models only require several tidal cycles. The reason for this is the relative speed of information propagation in both systems. In the hydrodynamic model, the information consists of small changes in sur face elevation that propagate along the estuary as elemental surges with speeds of u±/gf-f whereas in the mass transport model, the I information con sists of changes in concentration (finite or otherwise) that propagate at speeds of u (if dispersion is ignored). This is discussed further in Sec tion C.3. Thus, in applying a mass transport model to determine the quasi-steady state response of an estuary for some waste loading condition, it is necessary to run the model for a succession of tidal cycles. In solving the hydrodynamic and tidally varying mass transport equations for the Fra'ser River Estuary, the freshwater discharge and tidal conditions are assumed to be quasi-steady. Thus, the model estuary "sees" a constant freshwater inflow and a succession of identical tides. In actual fact, "slow" variations occur in the freshwater inflow and tidal 39 conditions, and the estuary may never achieve strict quasi-steady state conditions. However, as the residence time of the Fraser River Estuary is only four to six days, it seems reasonable to treat the freshwater and tidal conditions as constant for this period of time. 3.1.3 The Model Estuary of the Hydrodynamic Equations. The Fraser River Estuary is described in detail in Appendix B. The model estuary, as described byfthe hydrodynamic equations, differs to some ex tent from the real estuary because of the assumptions made in deriving the hydrodynamic equations. The model estuary consists of the same seven principles as the real estuary, although the Canoe Pass Area has been simpli fied into a single channel, as is shown in Figure B.4. The channels of the model estuary are assumed to be rectangular in cross-section (equal storage and advective widths), and the salinity is assumed to be everywhere - zero. Because the salt wedge has been ignored, the predicted velocities in the lower reaches of the estuary are probably somewhat low. The freshwater dis charge in the model estuary is constant, and the tidal rise and fall of the water surface is quasi-steady. 3.2 THE TIDALLY VARYING MASS TRANSPORT MODEL 3.2.1 Method of Solution. The tidally varying mass transport equation was solved by a characteristic finite difference method. In this method, the equation is solved along the characteristic curves of the advec tive transport process, as is discussed in Appendix C. The accuracy and directness of solution were instrumental in this choice; there is no numeri cal dispersion (see Appendix C) and the advection process is simulated 40 directly and independently of the dispersive process. It was not possible to use the usual technique's of dye patch studies or measurement of existing tracers to verify the mass transport models. The large dilutional capacity of the river made extensive dye patch studies too expensive, and as the seven channels of the estuary are essentially unpolluted [Benedict, et al. } 1973], there are no suitable existing pollutional tracers. Because of the high degree of stratifica tion and the complicated nature dynamics of the salinity intrusion (see Section B.2.4), salt could not be used as a tracer. Thus, it was not possible to verify the mass transport models by a direct comparison of predicted and observed concentrations. An attempt was made to "verify" the underlying mass transport processes of advection and dispersion inde pendently of each other, as is discussed in Section 4.2. The method of characteristics leads to a natural separation of the advective and disper sive transport processes, and this was one of the reasons for using this method to solve the tidally varying mass transport equation. The one-dimensional mass transport equation is derived in Appendix A and the assumptions in its derivation are discussed there. The equation was transformed into Lagrangian or characteristic form (see Section A.6) to give it. =' - n — = r"5 [EAr—} + 2 S (3>5) dt Adx 3x 1 i=l Equation (3.5) was then separated into its component dispersive and source-sink effects to give 41 dc dt 1 3 A 3x (3.6) and dc_ dt n Z S. (3.7) 1=1 1 Initially, a grid of moving points was assigned throughout the estuary. In the advection step, a finite difference form of Equation (3.4) was used to advect the points along the estuary for a time increment, the or through converging junctions. In the dispersion step, an explicit finite difference form of Equation (3.6) was used to adjust the concentration of the moving points for the effects of dispersion during the time increment. Finally, in the source-sink step, Equation (3.7) was used to adjust the concentration of the moving points for the effects of the source-sink pro cesses during the time increment. (The source-sink step can be done analy tically — for a BOD-DO system, Equation (3.7) is the usual Streeter-Phelps equation). The sequence of these three steps was then repeated for the next time increment, and so the solution progresses through time. At the boundaries of the estuary (the sea, Pitt Lake and Chilliwack), moving particles were added to and removed from the estuary as dictated by the advective boundary conditions. A time increment of one hour was used in solving the tidally varying mass trans port equation, and at the end of each hour,the concentrations were extrapolated off the grid of moving points onto the 5,000 foot fixed grid of Figure 3.1. The method of solution is very similar to the combined method of Gardiner et al. concentration of the points being adjusted as they passed effluent outfalls 42 [1964], except that whereas Gardiner et al, extrapolate the moving points onto the fixed grid, the moving points in the Fraser River Estuary solution always remain on their respective characteristics and only their concentra tions are extrapolated onto the fixed grid. Details of the solution of Equations (3.4), (3.6) and (3.7) are given in Appendix D. An advantage of the characteristic method of solution is that additional moving points can be added to the estuary to more closely de fine regions of rapid variation in concentration, as occur at times of slackwater at effluent outfalls, and where the variation is slow unneces sary moving points can be removed. Once the moving points are in the estuary, their subsequent positions are determined by the advective transport process, and they are not uniformly spaced along the estuary in'the dispersive step. Both an explicit and implicit finite difference scheme were investigated for the dispersive step. The explicit scheme was slightly slower, but was chosen for reasons of simplicity (see Appendix D). The stability requirements of the explicit scheme are discussed in Appendix C, and generally the time in crement for stability was less than one hour; the dispersive step then con sisted of a number of "internal" iterations within the basic time increment of one hour, as discussed in Appendix D. 3.2.2 The Model Estuary. For the sake of convenience, the model estuary of both the tidally averaged and tidally varying mass transport equa tions has been simplified to the three major channels of the real estuary, the Main Arm - Main Stem, the North Arm and the Pitt System, as shown in Figure 3.2. The stations along these channels are the same as those of the Pitt Lake '1/ fir _ Main Stem Figure 3.2 Simplified Model Estuary of the Mass Transport Equations 44 hydrodynamic model. It is possible to include all seven channels of the real estuary, but it would have made the computer programming of the tidally varying mass transport equation considerably longer and more involved. Because of their short length and smaller flows, the concentration profile along the minor channels is essentially determined by the profile along the major channels, and the approximation of ignoring the minor channels seems reasonable. In the model estuary, the dispersion process is assumed to be Fickian in its entirety and the cross-sectional mixing is assumed to be complete at the junctions. The salinity of the model estuary is zero, but because of the salinity intrusion in the real estuary, the predicted velocity field and advective transport in this region is probably under estimated. The estuary is highly stratified in the region of the salinity intrusion and little mixing occurs between the fresh and ^.saltwater. This will tend to minimize any chemical or biological effects the saltwater may have on the dissolved substance in question. As in the model estuary of the hydrodynamic equations, the freshwater discharge is constant and the tidal conditions are quasi-steady. 3.3 TIDALLY AVERAGED MASS TRANSPORT MODEL 3.3.1 Method of Solution. The one-dimensional, tidally averaged mass transport model of Thomann [1963] was used to determine the steady state tidally averaged response of the Fraser River Estuary for various waste loading conditions. (As the model is tidally averaged, the response is steady-state rather than quasi-steady state). In this model, the estuary is divided into 45 a number of longitudinal segments or "boxes", as shown in Figure D.2. Each segment is assumed to be completely mixed, and the tidally averaged advection and dispersion processes transport dissolved substance into and out of each segment. If the tidally averaged transport processes and efflu ent discharges remain steady, a mass balance about segment i gives a. . .-c + a. . c. + a. ...c' = W . (3.8) i,a-l i-l 1,1 i i,i+l l+l i where and c^ is the concentration in segment i; is the mass of effluent discharged into segment i per tidal cycle; a. . ^ is a coefficient accounting for the tidally averaged ' transport of substance between segments i and i-l, similarly for a^ ^+^; a. . accounts for dispersion out of segment i and any sink ' effects the substance undergoes. Details of this eguation and its coefficients are given in Appendix D. An equation similar to (3.8) can be written for each of the n segments of the estuary, and in matrix notation the system of equations can be written A • C = W (3.9) where | is a (n x 1) column matrix of the tidally averaged waste loads into.each segment per tidal-cycle; C is a (n x 1) column matrix of the tidally averaged concentration in each segment; 46 and ^ is a (n x n) tri-diagonal matrix containing the tidally averaged transport terms and any sink effects. Thus, the steady state tidally averaged response of the estuary is described by a system of n simultaneous linear equations. Equation (3.9) represents the concentrations of a conservative substance or a substance-^undergoing first order decay, such as BOD. To investigate the steady-state BOD-DO res ponse of an estuary, a system of equations similar to (3.9) and coupled to it are obtained for the DO concentrations (see Thomann [1971] for details). Thomann's steady-state solution to the one-dimensional tidally averaged mass transport equation is essentially a fixed grid finite differ ence solution similar to that of an implicit difference scheme (see Appendix D). There is no stability requirement for the steady-state solution (see Section C.3), but there is a non-negativity requirement that imposes rela tive limits on the sizes of the advective and dispersive transport processes. This is discussed in Section D.3, and if this requirement is violated, the concentration in a segment becomes negative. Thomann's steady-state -solution does not suffer from numerical dispersion, but for a rather unusual reason. Because of the fixed grid nature of his finite difference scheme, the advec tive process will not be correctly simulated and numerical dispersion should occur. However, from Section C.2, it is seen that for numerical dispersion to occur, it is necessary for the concentration at a fixed grid point to change with time. Thomann's steady-state solution admits no temporal changes, and thus no numerical dispersion occurs. Numerical dispersion does occur in his transient solution where the concentration at the fixed grid points does change with time, as Thomann [1971] recognizes. 47 3.3.2 The Model Estuary. The model estuary of the tidally aver aged equation consists of the same three major channels and stations as the model estuary of the tidally varying equation and cross-sectional mixing is assumed to be complete at the junctions. The tidally averaged model estuary has no%tides and has higher dispersion than its tidally varying counterpart. Thomann [1971] lists typical values of the tidally averaged dispersion coeffi cient. They range from 1-20 square miles per day with a mean value of about 10 square miles per day. These values seem high compared to the tidally varying values of dispersion coefficients, and apparently reflect the influence of tidal advection in determining the tidally average disper sion coefficient. As in the tidally varying solution, the influence of the salinity intrusion is ignored. CHAPTER 4 VERIFICATION OF THE HYDRODYNAMIC AND MASS TRANSPORT MODELS 4.1 THE HYDRODYNAMIC MODEL The verification of the hydrodynamic model is to ensure that the model estuary, as represented by the hydrodynamic equations, adequate ly reproduces the variation in water surface elevations and discharges (or advective velocities) observed in the real estuary. The "verification" consists of adjusting the friction factors and the cross-sectional widths and depths of the model estuary until an adequate fit is obtained between predicted and observed results. (In deriving the hydrodynamic eguations, the advective flow was assumed to be uniformly distributed over an assumed rectangular cross-section. In real situations, the cross-section is not rectangular and the flow is concentrated in the deeper sections. To compensate for these effects, it is necessary to adjust the cross-sectional geometry of the model estuary). 4.1.1 Data Available. The network of permanent tide-gauging stations throughout the estuary is shown in Figure B.ll, and provides a limited but adequate record of the water surface elevations of the estuary throughout the tidal cycle. Field measurements of velocities and discharges are sporadic and inadequate for verification purposes. The only existing data adequate for a complete verification of the hydrodynamic model under high tide-low flow conditions is due to Baines [1952]. In this study, water surface elevations were recorded at half-hour intervals at 43 stations 48 49 throughout the estuary for the high tide-low flow conditions of January 24, 1952. The freshwater discharge at Chilliwack was 36,500 cubic feet per second and the tidal range at Steveston was 11 feet. Baines used the method of cubature to estimate the tidally varying discharges at the 43 stations. 4.1.2 Water Surface Elevations for Low Flows. The high tide-low flow conditions of January 15, 1964 were used in an initial attempt to reproduce the recorded water surface elevations in the real estuary. The freshwater discharge at Chilliwack was 53,500 cubic feet per second and the tidal range at Steveston was 10 feet. The water surface elevations of the model estuary were found to be relatively insensitive to the effects of friction and cross-sectional geometry. In fact, the response of narrow sections with high friction was found to be equivalent to that of broader sections with low friction. This is illustrated in Figure 4.1 which shows the predicted responses for both types of section at Station No. 18 (New Westminster) on the Main Arm - Main Stem. The average values of widths, depths and Manning's "n" along the Main Arm —. Main Stem are also shown. Note that while Chezy's formula for friction was used in solving the hydro-dynamic equations, friction coefficients are reported as values of Manning's "n" (the values of n are for feet-second units). 4.1.3 Water Surface Elevations for High Flows. To separate out the independent effects of friction and cross-sectional geometry, an attempt was made to reproduce the water surface elevations of the high tide-high flow conditions of June 16, 1964. The freshwater discharge at Chilliwack was 463,000 cubic feet per second and the tidal range at Steveston was eight feet. It was thought that the greater velocities under high discharge 50 Observed ° o Predicted No.I (Narrow section,high friction) * x Predicted No.2 ( Broad section,low friction) CJ> O CO o o Q = 52,400cfs (Chilliwack ) Tidal Range of Steveston =10 , J Jan.15,19640^ 8 12 16 Hours 20 M Average values along Main Arm - Main Stem Results Average Mannings •i II n Average Width (feet ) Average Depth (feet ) Predicted No. 1 0.042 1350 30.6 Predicted No. 2 0.030 1870 22.9 Figure 4.1 Equivalent Stage Response for Two Different Types of Cross-Section 51 conditions would be more sensitive to frictional effects. This was found to be the case, the tidally averaged or mean water levels being very sen sitive to friction and essentially independent of width. However, the width was found to govern the range of the water surface fluctuations about the mean water level. The predicted and observed water surface elevations are shown in Figures 4.2 and 4.3, and the reproduction of the observed re sults is seen to be satisfactory. The gross cross-sectional values of depths and areas shown in Figures B.6 to B.8 were used in obtaining these results. The width was determined by dividing the area by the depth, and was a satisfactory compromise between the narrower advective sections and the broader sections more representative of the storage width. Manning's "n" varied from 0.022 in the lower reaches to 0.027 in the upper reaches. These values are reasonable, and indicate that the Main Arm - Main Stem is hydraulically smooth [Chow, 1959]. Under these high discharge conditions, the flow in the upper reaches of the Main Stem is steady, and the predicted water surface slope can be used to check that the hydrodynamic model satisfactorily reproduces the frictional effects. The predicted value of Manning's "n", as determined by the average area, depth and predicted water surface slope between stations 40 and 60 on the Main Stem, was 0.0269. This agrees closely with the actual value of 0.027 used in the hydrodynamic model over this section of the estuary. The "high flow" friction coefficients and cross-sectional geo metries were then used in a second attempt to reproduce the water surface elevations for the low flow conditions of January 15, 1964. To obtain a sat isfactory fit between the predicted and observed results it was necessary to I I I I 1 I I M 4 8 12 16 20 M Hours Figure 4.2 Observed and Predicted Stages for June 16, 1964 Q= 463,000 cfs (Chilliwack) Tidal range at Steveston s 8' June 16, 1964 20 8 I 6 Main Arm - Station No.22 ( Port Mann) 24 |~ Main Arm - Station No. 32 ( Port Hammond) 22 20 0 bserved o o Predicted 22 I" Pitt - Station No. 144 ( Port Coquitlam) OA o O O ° O ° rt o r> ~ 20<K°. • -2— ° n n n n n n n IL I 8 22 r Pitt - Station No. 162 ( Al vi n ) 20 I 8 o o M 8 12 Hours 16 20 M Figure 4.3 Observed and Predicted Stages for June 16, 1964 54 increase Manning's "n" by approximately 0.005 throughout the estuary. Under these conditions the value of Manning's "n" along the Main Arm-Main Stem varied from 0.027 in the lower reaches to 0.03 2 in the upper reaches, and once again these values are reasonable. The increase in fric tion that low flows experience relative to higher flows also seems reason able. During the falling limb of the annual hydrograph sediment: is deposited in the lower channels of the estuary [Pretious, 1972], and thus the low winter flows experience bed forms defined by the higher summer flows. It is expected that the relative roughness of these bed forms will be greater for the low flows — and hence the greater friction. 4.1.4 Verification for the Conditions of January 24, 1952. Having made an independent check on the friction coefficients and effective cross-sectional geometries, an attempt was made to reproduce the recorded water surface elevations and cubature discharges of Baines [1952]. With minor adjustments to the friction coefficients, the water surface elevations and discharges shown in Figures 4.4 to 4.11 were obtained. The predicted dis charge curves were found to be insensitive to changes in the cross-sectional geometries, but relatively sensitive to changes in the friction coefficients. In certain sections of the estuary, the low-flow friction coefficients of the section have been slightly lowered to obtain a better fit between the pre dicted and observed peak discharges. Generally, the fit between the predicted and observed water surface elevations is satisfactory. In the upper reaches of the Main Stem, the range of predicted water surface elevations is somewhat greater than that observed, and is probably due to "ending" the model estuary before the tidal rise and 55 Q = 36,500 cfs (Chilliwack) M 4 8 12 16 20 M Hours Figure 4.4 Tidally Varying Stage and Discharge - Main Arm Q = 36,500 cfs (Chilliwack) m-jTr, x.r.i. x. January 24, 1952 Tidal Range at Steveston =11' J ' 56 E o •D U> O <» •o c o CO OJ a> 'T 12 I I 10 9 8 7 6 5 4 3 2 Main Arm - Stat.No.14 ( St. Mungo Cannery) O 3 O o 240 -200 -160 -120 -80 -40 0 40 80 120 160 </) 200 ° O O 240 O a> a> 14 cn ° 13 CO I2h I I 101 9 8 7 6 5 4 3 \ // / o Z3 O o _l Lu z Main Arm - Stat. No.18 (New Westminster R.R.Bridge) I I I I 240 £ in -200 a -160 -120 -80 -40 0 40 80 120 180 200 M 8 12 Hours 16 20 M Figure 4.5 Tidally Varying Stage and Discharge - Main Arm Q = 36,500 cfs (Chilliwack •- q> cn I i-o M 4 8 12 16 20 M Hours Figure 4.6 Tidally Varying Stage and Discharge - Main Arm e O T3 W TJ O c o CO CO q> cu o> o CO 15 14 13 12 I 1 10 9 8 7 6 5 4 3 2 Q = 36,500 cfs (Chilliwack) m-j-ir, January 24, 1952 Tidal Range at Steveston =11' J ' Observed Stage Predicted Stage Cubature discharge IBaines) Predicted discharge Main Arm - Stat. No.40 ( Whonock) J_ 20 M 8 12 16 Hours Figure 4.7 Tidally Varying Stage and Discharge - Main Arm * "I o -i U-Z 58 140 H-120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 u O O O I -14,0 o> t_ o -120 x: o «A -100 Q -80 - 60 -40 -20 0 20 40 60 80 1 0 0 Q = 36,500 cfs (Chilliwack) 59 m-j-ir, j-ox. n.} January 24, 1952 Tidal Range at Steveston =11' 11 M 4 8 12 16 20 M Hours Figure 4.8 Tidally Varying Stage and Discharge - North Arm 15 't 13 h 12 I I I— 10 Q = 36,500 cfs (Chilliwack) Tidal Range at Steveston = 11' } January 24, 1952 •'/ Cubature discharge t Baines ) \ 1/ Predicted discharge \\J Middle Arm - Stat. No. 133 ( Airport Floats ) 60 •30 25 20 - 15 - 10 5 I. 10 -J 15 20 H25 30 35 Canoe Pass - Stat. No. 124 ( Canoe Pass) o -I u. I-o 5 o _i • b. Z -70 -60 -50 -40 - 30 - 20 - 10 0 10 20 30 40 50 M 4 8 12 16 20 M Hours Figure 4.9 Tidally Varying Stage and Discharge - Middle Arm, Canoe Pass Q = 36,500 cfs (Chilliwack) Tidal Range at Steveston = 11' } January 24, 1952 2 12 I I 10 9 8 7 6 5 4 3 2' 15 14 h-13 12 11 10 9 8 7 6 5 4 3 \\ Pitt - Stat.No. 144 (CPR Bridge Coquitlam) 61 -120 -100 -80 -60 - 40 - 20 0 20 40 60 80 100 120 • j Pitt - Stat. No. 151 ' • (Gilley's Quarrey) - 140 H-120 -100 - 80 - 60 - 40 - 20 0 20 40 60 80 100 u O O O Urn a o M M 8 12 16 Hours 20 M Figure 4.10 Tidally Varying Stage and Discharge - Pitt River Q = 36,500cfs (Chilliwack ) Tidal range at Steveston = ll' Observed Stage Predicted Stage M 4 8 12 16 20 M Hours Figure 4.11 Tidally Varying Stage — Pitt Lake j Jan. 24,1952 63 fall of the water water surface has dissipated. The fit between the pre dicted discharges and the cubature discharges of Baines is not as good as the fit between predicted and observed water surface elevations. There seems to be a phase difference between the discharges, the predicted dis charge tending to occur earlier than the cubature discharge. This phase difference was insensitive to changes in cross-sectional geometry and fric tion coefficients. The cubature discharge curves behave somewhat erratically between the times of high-high-water and low-high-water. It is not certain whether this behaviour represents the true discharge response of the estuary, or is due to inaccuracies in the cubature discharge calculations. The tidally averaged discharges were found to satisfy continuity at the Main Arm - North Arm junction, the North Arm - Middle Arm junction and the Main Arm - Canoe Pass junction. A tidally averaged mass balance was made over the entire estuary and was found satisfactory. The results of this mass balance are shown in Table 4.1. To sum up, the behaviour of an estuary is determined by the inter action of the effects of river flow, tides, bed forms and cross-sectional geo metry. The freshwater discharge of the Fraser River Estuary undergoes a large annual fluctuation (see Appendix B) and consequently the bed forms of the var ious channels, and possibly their cross-sectional geometries, are in a contin ual state of dynamic equilibrium with the freshwater flows and tides. Bearing this in mind, and the fact that Baines obtained his results some 20 years ago, and additional dykes, training walls, etc. have been constructed in the interim, it is impossible for the simple one-dimensional hydrodynamic model to exactly reproduce the behavior he observed in the estuary. The 64 TABLE 4.1 TIDALLY AVERAGED MASS BALANCE OF PREDICTED DISCHARGES FOR JANUARY 24, 1952 CHANNEL PREDICTED VALUES DISCHARGE IN (1000*s cfs) MAIN STEM 36.5 DISCHARGE OUT (1000*s cfs) MAIN ARM NORTH ARM MIDDLE ARM CANOE PASS 27.8 2.9 0.5 1.2 TOTAL 32.4 INTERNAL STORAGE (1000's cfs) PITT SYSTEM OTHER 2.9 M).0 TOTAL 2.9 results shown in Figures 4.4 to 4.11 represent the "best over-all fit" between predicted and observed water surface elevations and discharges. As such, the predicted results are close to the limit of resolution of the one-dimensional hydrodynamic model, or in other words, the predicted results are the "best" the model is capable of. 65 4.2 THE TIDALLY VARYING MASS TRANSPORT MODEL Because of the lack of field data (see Section 3.2.1), a thorough verification of the tidally varying mass transport'-, model was not possible. However, the model is computationally stable, free from the effects of numerical dispersion and sufficiently flexible to be adjusted to field data when available. It reproduces the advective transport process correctly and demonstrates good agreement with the standard analytical result for the dispersion of a slug load. To obtain preliminary notions of the tidally varying behaviour of the estuary, other peoples' results have been used to estimate the coefficients of longitudinal dispersion. 4.2.1 Advective Transport. In the tidally varying model, the advective transport was simulated by a grid of moving points on the advec tive characteristics, as described in Chapter 3. This process directly simu lates the advective transport occurring in the actual estuary and does not suffer from numerical dispersion (see Appendix C). The hydrodynamic model was used to predict the velocities throughout the estuary at half-hour inter vals and these velocities were then used to advect the moving points along. Thus, the tidally varying mass transport:: model will simulate the advective transport as accurately as the hydrodynamic model simulatesv the velocities, and the verification of the hydrodynamic mode] can also be regarded as a partial verification of the advective transport process. The predicted ad vective transport was found to be relatively insensitive to inaccuracies as in the advective velocities. This is illustrated in Figure 4.12 which shows the advection of a slug load down the Main Arm - Main Stem of two different model estuaries. The freshwater discharge and quasi-steady tidal conditions are the same for both estuaries, but the higher friction of one estuary January 24, 1952 - Flow at Chilliwack « 36,500cfs 20 h Predicted No.I ( (Marrow Section with high friction) S— Predicted No. 2 ( Broad Section with low friction) AVERAGE VALUES ALONG MAIN ARM - MAIN STEM RESULTS AVERAGE MANNINGS "N" AVERAGE WIDTH ( feet ) AVERAGE DEPTH ( feet) PREDICTED No. 1 0.042 1 3 50 30.6 PREDICTED No.2 0.0 30 18 70 22. 9 20 40 60 Hours 80 100 120 Figure 4.12 Advection Of A Slug Load Down The Main Arm - Main Stem CTl CT* 67 distorts its velocity field relative to the velocity field of the other. The only difference between both sets of results is in the tidal excursion in the lower reaches of the estuaries, and this is not that significant. When the tidally varying model was run with zero dispersion, the dilution and advection of effluent through the junctions was found to be correct (it being recalled that cross-sectional mixing is assumed to be complete at the junctions). The discrepancy in the mass balance? over any tidal cycle was found to be less than 6% and was due to the discrete spacing between the moving points. 4.2.2 Dispersive Transport. A slug load was introduced at Station No. 50 on the Main Stem to verify the capability of the tidally varying mass transport model to simulate the dispersion process. The fresh water discharge at Chilliwack was 36,500 cubic feet per second and the dis persion coefficient was set equal to 500 square feet per second and assumed constant in x and t. The analytical solution for the dispersion of a slug load is given by [Fischer, 1966a] c = ,m exp{- —} Vwvt 4Et where M is the mass per unit area introduced into the flow. For the given conditions, this reduces to 136 r 3.5X2-, c = "Tt exP{ —} where c is the concentration in milligrams per litre; 68 t is the elapsed time in hours; and x is the distance either side of the mean value in station co-ordinates (5,000 foot segments). Figure 4.13 shows the form of the slug inputs for the analytic and tidally varying solutions. (The effluent discharge in the tidally varying model sim ulates continuous discharges, hence the form of the slug input in Figure 4.13). The predicted and analytic solutions for the dispersion of the slug load are shown in Figure 4.14 and the agreement between both sets of results is good. The higher peak concentrations of the analytic solution are due to the higher initial concentrations of the analytic slug load. It is recalled that the dispersion process is assumed to be Gaussian in its entirety, and this is illustrated by the results of Figure 4.14. The estimation of the coefficients of longitudinal dispersion is discussed in detail in Appendix F. In obtaining the tidally varying results of this study, the coefficient of longitudinal dispersion is assumed to be given by E = (6 + S£)y E = ay U 0 < t < T t > T (4.1) and where = 0.06u (4.2) E/ Yr u* an<^ u are the "instantaneous" values of the respective parameters during the tidal cycle; t is the time that has elapsed since a "parcel" of the effluent was discharged into the estuary; 2 T = b /e and is the time scale of lateral mixing (an edge dis charge is assumed); 69 Discharge at Station No. 50 Q = 36,500cfs (Chilliwack) Tidal Range at Steveston = 11 feet J \ Jan.24,1952 Mass = 1650x10 lbs. 150 Slug load of Analytic Soluti on Slug load of Predicted Solution 21 22 Time ( hours ) Figure 4.13 Slug Inputs for Analytic and Predicted Dispersion Solutions Q = 36,500 cfs (Chilliwack) ) jan 24 1952 Tidal Range at Steveston = II feet J ' 3 HOURS AFTER INJECTION Predicted + Analytic 136 . . f-3.47x" C = exp.{ —} t in hours x in station coordinates J L 50 Stations 55 40 r-20 0 35 6 HOURS AFTER INJECTION 45 Stations 50 18 HOURS AFTER INJECTION 45 Stations Figure 4.14 Dispersion Of A Slug Load in The Main Stem 71 and a is a measure of the effects of lateral velocity gradients on the dispersion process and is tabulated in Tables F.l. Equations (4.1) allow for the increasing contribution of the lateral velo city gradients on the dispersion process as the effluent spreads across the section and the variation of the coefficient during the tidal cycle. These aspects are discussed in Appendix F and it is seen that that the peak tidally varying concentration is quite sensitive to assumptions about the form and magnitude of dispersion coefficient. Although simplistic,the relationships of Equations (4.1) are thought to be a reasonable approximation of what actually occurs in the estuary. Because the tidally varying mass transport equation hastbeen solved along the advective characteristics, the position of each effluent "parcel" and the time that it has spent in the estuary is known. (This in formation is "masked" in a fixed grid solution). Consequently, the tidally varying mass transport model can account for time dependent behaviour of in dividual "parcels" of effluent, as is assumed in Equations (4.1). 4.3 THE TIDALLY AVERAGED MASS TRANSPORT MODEL The tidally averaged dispersion coefficient includes the effects of advective transport due to tidal flow reversal and longitudinal dispersion. As such, it has no real physical meaning and the verification of a tidally averaged mass transport model essentially consists of forcing the model to fit field data. According to Holley et al. [1970] there is no way of making an a -priori estimate of this coefficient. Ward, and Espey [1971] suggest using the results of a tidally varying model to estimate the coeffi-72 cient. In the absence of field data, this latter approach has been adopted to obtain estimates of the tidally averaged dispersion coefficients for the Fraser River Estuary. Thomann [1971] lists values of the tidally averaged dispersion coefficient for various estuaries. They range from 1-20 square miles per day with a mean value of around 10 square miles per day. Figure 4.15 shows the influence of the tidally averaged dispersion coefficient on the predicted concentrations for a steady, continuous discharge of a conservative effluent at Station No. 40 on the Main Stem of the Estuary. (The concentrations are standardized by dividing by the tidally averaged concentration based on the total mass of effluent discharged during the tidal cycles and the freshwater flow at Chilliwack. This is plotted as the y , the subscript signifying ta that the results have been obtained from the tidally averaged model). The only significant difference between the results is in the upstream excursion. The concentration gradients are small downstream of the effluent discharge, and the dispersion has little effect. It is recalled that the results of Figure 4.15 have been obtained for a conservative effluent. If the concen tration of the effluent decays with time, as with BOD, there will be concentra tion gradients downstream of the discharge point, and the downstream concentra tions will also depend on the values of the tidally averaged dispersion coeffi cient (although the downstream concentrations^are not very sensitive to the values of the dispersion coefficient). The decrease in concentrations around Stations Nos. 18 and 24 is due to the effects of the junctions and introduces an error of about 10% into the predicted concentrations. It is not known why this effect occurs. Q = 36,500cfs (Chilliwack) )Jan.24>|952 Tidal Range at Steveston =11 feet J l.2r-1.0 -0.8 -0.4 -0.2 0 E = lOsquare miles per day 10 20 30 Stations along Main Arm 40 50 Main Stem Figure 4.15 Influence Of The Tidally Averaged Dispersion Coefficient Oh Predicted Tidally Averaged Concentrations 74 Figures 4.16 and 4.17 show the maximum upstream excursion on the flood flow for effluent released along the various channels. In the lower reaches of the Main Arm it is assumed that the tidally averaged dispersion coefficient is equal to 20 square miles per day, and this value is reduced along the Main Arm - Main Stem according to the results of Figure 4.16. Along the North Arm the tidally averaged dispersion is assumed to equal 10 square miles per day, and along Pitt River it is assumed to equal 30 square miles per day (the rapid decrease in the upstream excursion along the Pitt River, as shown in Figure 4.17, is due to the drastic decrease in advective velocity in Pitt Lake). Although the absolute values of the tidally averaged dispersion may not be correct, and this is discussed in Chapter 5, the relative values along the estuary should be reasonable. Maximum Upstream Excursion (Stations) Q - 36,500cfs (Chilliwack) )Jan.24il952 Tidal Range at Steveston = 11 feet J Effluent moves into Pitt Lake 01 c o o CO o 2 8 c o « 6 i_ 3 O X LU E 4 t_ to Q. Z> E 2 3 E 0 Effluent moves through N.Arm - M.Arm Junction J L Pitt Lake "00 110 118 140 150 Position of Effluent Discharge along North Arm and Pitt River ( Station Numbers ) Figure 4.17 Maximum Upstream Excursion During Flow Reversal In The North Arm And Pitt River CHAPTER 5 COMPARISON AND DISCUSSION OF RESULTS The three models discussed in the last two chapters are now used to obtain'a.preliminary indication of the significance of tidal effects on predicted concentrations in the Fraser River Estuary. The flow and tidal conditions of January 24, 1952 were used in this investigation. These are typical high tide - low flow conditions and were used in the verification of the hydrodynamic model. The freshwater discharge at Chilliwack was 36,500 cubic feet per second and the tidal range at Steveston 11 feet. A conservative effluent and an edge discharge were assumed. The hydrodynamic model was used to obtain the tidally varying velocities and cross-sectional areas along the various channels of the estuary. These values, together with estimated dispersion coefficients, were then used to predict the tidally varying concentrations. The freshwater discharge at Chilliwack, the tidally averaged areas and estimated dispersion coefficients were used to pre dict the tidally averaged concentrations. Before presenting the results of this investigation, the assumptions of the mass transport models are briefly recalled: both models are one-dimensional and predict the cross-sectionally averaged concentrations; the cross-sectional mixing in both models is assumed to be complete at the junc tions; the influence of the saltwedge is ignored in both models; the dispersion process of the tidally varying model is assumed to be Gaussian in its entirety; and in the tidally averaged model, all effects of the tides are lumped into the tidally averaged dispersion coefficient. 77 78 5.1 THE EFFECTS OF LATERAL DISPERSION In the initial period before cross-sectional mixing is complete, there are significant lateral concentration gradients across the estuary, and the peak cross-sectional concentration will be considerably higher than the value predicted by either mass transport model. While this effect has not been quantitatively assessed, it is noted that the tidally varying mass transport model can easily be adapted to give an approximate estimate of this effect. The lateral concentration profile depends on the coefficient of lateral dispersion, the position of the effluent outfall in the estuary cross-section and the period of time that a "parcel" of effluent has spent in the estuary [Ward, 1972]. Because the tidally varying mass transport model has been solved along its advective characteristics, the position of each effluent "parcel" and the time that it has spent in the estuary is known, and the predicted concentrations can be adjusted to account for the effects of lateral dispersion. (This approach was used in allowing for the assumed time dependent variation in the coefficient of longitudinal dispersion). In the lower reaches of the estuary, there are three distinct periods in each double tidal cycle of 25 hours; a period of strong flood flows; a per iod of weak ebb and weak flood flows and a period of strong ebb flows. This is apparent in Figures 4.4 to 4.10 and the duration of each of the three per iods is seen to be approximately eight hours. The cross-sectional mixing is principally due to the effects of secondary flows, as discussed in Appendix E, and will be greatest during the periods of strong ebb and flood flows. If cross-sectional mixing is not complete at the Main Arm - North Arm junction during the strong ebb flow, or at the Main Stem - Pitt River junction during 79 the strong flood flow, the effluent will not be advected through these junc tions according to the simple flow balance of the one-dimensional mass trans port models. The interaction of the time dependent behaviour of the lateral mixing process with the advection of effluent through the junctions is very involved and can only be reliably determined from field studies. If suffi cient field data were available, it may be possible to empirically allow for this effect in the tidally varying model. 5.2 THE INITIAL DILUTION OF EFFLUENT When effluent is discharged into an estuary, the tidally varying flows cause a variation in its initial dilution during the tidal cycle. This process is described in Section 1.3 and is seen to generate concentration spikes at times of slackwater. Figure 5.1 shows typical variations in the initial dilution during the double tidal cycle at four stations along the estuary. During a double tide cycle there are four slackwaters in the lower reaches of the estuary (see Figures 4.4 to 4.10) and the resulting concentration peaks at Stations Nos. 10, 22 and 102 are apparent. At Station No. 50, the tidal flows are significantly smaller in relation to the freshwater flow, and their combined interaction results in only two effective peaks during the double tidal cycle. Because of the smaller flows along the North Arm, the tidally varying concentrations at Station No. 102 are significantly higher than the values along the Main Arm - Main Stem. The tidally averaged dilution is also shown in Figure 5.1, the higher values at Stations Nos. 10 and 102 reflecting the division of the freshwater discharge between the North Arm and Main Arm. Q= 36,500cfs ( Chilliwack ) \jan. 24,1952 8° Tidal Range at Steveston = llfeetj 4 2 0 4r c o 0 Station No. 10 12 18 24 (Hours) 12 - Station No. 22 —•—*~v 1 ' i ^ r— 18 24(Hours) Station No.50 18 24(Hours) Station No.102 6 12 Figure 5.1 Initial Dilution At Station Nos. 10, 22, 50 and 102 24(Hours) 81 The peak tidally varying concentrations in Figure 5.1 are seen to be about twice the tidally averaged values. The minimum tidally varying concentration occurs during the strong ebb and flood flows, and is reflected as the flat sections of the curves of Figure 5.1. At Stations Nos. 10 and 22, this minimum concentra tion is approximately half the minimum value at Station No. 50. This is due to the presence of Pitt Lake; discharges downstream from the Main Stem-Pitt River junction are diluted by the combined effects of the tidal prisms of the Pitt system and the Main Stem above the junction, whereas discharges upstream of the junction are diluted only by the latter. 5.3 PEAK EFFLUENT CONCENTRATIONS As well as the effects of variation in initial dilution, certain slugs of water are dosed with effluent several times due to the effects of tidal flow reversal. This process is described in Section 1.3 and also gen erates concentration spikes. However, the effects of variation in initial dilution and multiple dosing combine to produce a compound spike whose concen tration is higher than that of the individual component spikes. This is well illustrated by an effluent discharge at Station No. 50. The variation of velocity at Station No. 50 is shown in Figure E.2; note the slackwaters at hours 4 and 5 and that the minimal velocity at hour 6. This is the reason for the increased concentration at hours 4 and 5 and the spike at hour 6 in Figure 5.1. Now consider the behaviour of the water mass around Station No. 50 during this time. During hour 4 a slug of water moves downstream past the outfall and is dosed, during hour 5 it moves back upstream and is dosed again 82 and finally during hour six it moves downstream past the outfall is dosed yet again. The result of this is that the concentration in this slug of water is 12 times greater than the tidally averaged value as is shown in Figure F.2 for the curve with zero longitudinal dispersion. This spike is composed of the three subspikes generated at hours 4, 5 and 6, their indivi dual magnitudes being 3, 3 and 6 as shown in Figure 5.1. The spikes due to the effects of variation in initial dilution and multiple ?dosing combine because both effects occur around the same phase of the tide. A particularly sensitive portion of the double tidal cycle is the period of weak ebb and flood flows where there are three slackwaters in the space of eight hours. The velocities between these times of slackwater are relatively low (see Figures 4.4 to 4.10). Once a spike is generated by the effects of variation in initial dilution, multiple dosing or a combination of both, its peak concentration begins to be reduced by the longitudinal dispersion process. This is illus trated in Figure 5.2 for a steady effluent discharge at Station No. 50. The compound spike is generated at hour 6 with a peak concentration that is 11 times greater than the tidally averaged value. (The decrease from the pre vious value of 12 is due to the effects of dispersion on the subspikes gen erated at hours 4 and 5). However, as the spike is advected down the estuary, the dispersion process reduces the peak concentration as shown. The results of Figure 5.2 seem to indicate that mass is not conserved since the area under the spike is not constant. This is only an apparent effect due to extrapolating the concentrations off the advective characteristics onto the standard 5,000 foot space grid, as described in Section 3.2. In actual fact, the initial base width of the spike is only 500-800 feet, which is too fine to be resolved 12 Q = 36,500cfs at Chilliwack ^ , oyl lftC.0 J ^ . I Jan. 24,1952 Tidal Range at Steveston II J 10 8 Vtv 6 o1— Stations along Main Stem Figure 5.2 Dispersion Of A Concentration Spike 1 hour 03 CO 84 by the fixed grid. (Note that the spike is correctly resolved on the char acteristics) . The most significant effects of variation in initial dilution and multiple dosing occur in the first 1-2 tidal cycles after a "parcel" of effluent has been discharged into the estuary. This is due to the high initial concentration gradients of the spike and assumed time-dependent increase in the coefficient of longitudinal dispersion (see Appendix F). These aspects are apparent from Figure 5.2 and are also illustrated in Appendix F. For the steady discharge of a conservative effluent, the peak tidally varying concentration was found to be from 1 to 10 times greater than the value predicted by the tidally averaged model. These results are summarized in Table 5.1" and were obtained for a single effluent discharge at the designated station. (The peak tidally varying discharge occurs at the discharge point). The results of Table 5.1 are sensitive to the magni tude arid assumed temporal variation of the coefficient of longitudinal dis persion, as discussed and illustrated in Appendix F. In obtaining the results of Table F.l, the coefficient of longitudinal dispersion was assumed to be as given in Equation (4.1). The peak cross-sectional concentration will be significantly higher than these values. 5.4 UPSTREAM EFFLUENT TRANSPORT In a tidally varying mass transport model, effluent is transported upstream by the effects of dispersion and advection during tidal flow reversal. Because of the large tidal flows, advective transport is the more important CHANNEL Main Arm - Main Stem North Arm Pitt River STATION 8 14 22 30 40 50 60 106 116 142 150 ( Ttv/ *Xa) max. 5.4 3.8 3.8 10.7 5.0 10.9 1.0 7.2 8.2 1.2 1.2 Table 5.1 Ratio of Peak Tidally Varying Concentration to Tidally Averaged Value At Point Of Effluent Discharge 86 component for the Fraser River Estuary, as is apparent from Figures 4.16 and 4.17. In a tidally averaged model, upstream transport is by the tidally averaged dispersion process. Figures 5.3 and 5.4 show the predicted tidally varying and tidally averaged concentrations at Stations upstream from effluent discharges. In all cases, the tidally averaged value is significantly less than the tidally varying values. A better fit between the results could be obtained by increasing the values of the tidally averagedMispersion coeffi cient. The results of Figure 5.4 are for simultaneous effluent discharges at Stations Nos. 10 and 14. Once again, the effects of variation in initial dilution and multiple dosing cause a concentration spike at both discharge points, spike A being generated at Station No. 10 and spike B at Station No. 14. Nbte that between hours 6 and 17 the tidally varying concentration at Station No. 6 is zero due to uncontaminated seawater moving upstream on the flood tide. The concentration boundary condition at the sea was assumed to be zero, a more realistic boundary condition being possibly as shown. The tidally averaged concentration is not a good indication of the tidally varying response upstream of Station No. 14, between the two stations or downstream of Station No. 10. 5.5 CHANNEL INTERACTIONS Consider now the effect of an effluent discharge in one channel on the water quality in another channel. The significance of this effect will depend on the proximity of the effluent discharge to the junction. If the discharge point is sufficiently far upstream from a junction, the effluent Discharge at Stat. No. 14 ( Ma in Arm) 6 12 18 24(Hours) Discharge at Stat. No.22 ( Main Arm) 6 12 18 24(Hours) Figure 5.3 Predicted Concentrations in One Channel Caused By Effluent Discharge In Another Channel Q = 36,500cfs (Chilliwack) Tidal Range at Steveston = Ufeet j Jan.24,1952 Station No.17 2r 24 (Hours) 0 8 r Station No.12 ' 1^ spi ke B - /*» 1 1 1 1 0 6 12 Station No.6 18 24( Hours) More Realistic Boundary Condition -Assumed Boundary Condition . spike A spike B 12 18 Figure 5.4 24( Hours) Predicted Concentrations Due to Two Effluent Discharges In The Main Stem 89 will be dispersed both laterally and longitudinally by the time it reaches the junction. Under these conditions, the tidally averaged concentration in both channels downstream of the junction will be a good measure of the tidally varying values. Of more interest is the case where the effluent discharge is downstream from a junction. This is illustrated in Figure 5.3 for interactions between the Main Arm and North Arm and between the Main Stem and Pitt River. Consider the effluent discharged at Station No. 22: at hour 8 a spike is advected up Pitt River past Station No. 144 on the strong flood flow and returns past Station No. 144 at hour 23 on the strong ebb tide, its concentration being reduced by the effects of dispersion dur ing its residence in Pitt River. For the two cases of Figure 5.3, upstream advection causes a slug of contaminated water to be fed into the other channel during each double tide cycle of 25 hours. The zero concentration in the North Arm at hour 22 marks the end of one slug and the beginning of another. In both cases, the tidally averaged concentration is not a good indication of the tidally varying values, although as remarked in a previous section, the fit can be improved by increasing the tidally averaged dispersion coefficient. 5.6 SUMMARY In this chapter the results from the tidally varying and tidally averaged mass transport models have been compared to obtain a preliminary indication of the significance of tidal effects on the predicted concentra tions in the Fraser River Estuary. 90 The tidally varying flows cause a variation in the initial dilution of a discharged effluent. This, and the effects of tidal flow reversal, independently generate concentration spikes. Because both effects occur around the same phase of the tide, the period of weak flood and ebb flows in each double tidal cycle being especially sensitive, the individual spikes combine-to generate a spike whose concentration was estimated to be from 1 to 10 times greater than the predicted tidally averaged concentration. Because of incomplete lateral mixing, the peak cross-sectional concentration will be considerably higher than these values. The effects of concentration spikes are most significant in the first 1-2 tidal cycles after their gener ation. After this period of time, the longitudinal dispersion is due to the effects of lateral concentration gradients and the concentration of the spike is rapidly reduced. The movement of water through the Fraser River Estuary is a com plex process that is influenced by the large tidal effects, the significant freshwater discharge (even at low flows), the various channels and junctions of the estuary and the presence of Pitt Lake. If an effluent discharge is sufficiently far upstream from the junctions, the spikes will be flattened and the effluent dispersed over the cross-section by the time it is advected to the junction. Under these conditions, the predicted tidally averaged be haviour at and downstream of the junction is a good estimate of the tidally varying behaviour, and both should be a reasonable approximation to the actual behaviour in the estuary. Effluent discharges in the Main Stem above Station No. 45 are expected to behave in this manner. Throughout the lower reaches of the estuary,,;including Pitt River, the predicted tidally averaged behaviour 91 is not a good indication of the tidally varying behaviour. Because of the effects of incomplete lateral mixing and junctions, as discussed in Section 5.1, the actual response of the estuary may be somewhat different from the predicted tidally varying behaviour. By adjusting the tidally averaged dispersion coefficients, it is possible to obtain a "best-fit" between the predicted tidally varying and tidally averaged concentrations. However, since the tidally averaged dis persion coefficient has no real physical meaning, there is probably no unique set of values for the estuary. Rather, the best fit dispersion co efficients will vary with the position and number of effluent discharges. CHAPTER 6 SUMMARY AND CONCLUSIONS Mass transport models are commonly used to investigate situations of existing or potential water quality conflicts in estuaries. Such models are used to predict the concentration of the offending substance through out the estuary and are essentially of two types: those that correctly allow for tidal effects and those that do not. While tidally varying models correctly allow for the effects of the tides, their development and application involves significantly more work than for their tidally averaged counterparts. In tidally averaged models, all tidal effects are lumped into the tidally averaged dispersion coefficient, and while such models are relatively easy to develop and apply, they give no indication of the significance of tidal effects on the predicted concentrations. In this study, both a tidally averaged and a tidally varying mass transport model have been developed and applied to the Fraser River Estuary. The principal object of this study was to assess the ability of the tidally averaged model to describe the tidally varying concentrations. This was investigated by comparing the predicted concentrations from both models for assumed effluent discharges. The tidally averaged model used in this study is due to Thomann [1963]. The tidally varying model was developed from first principles, and during its development it was necessary to consider the problems of numerical dispersion, stability, the significance of lateral 92 93 dispersion and the time dependent behaviour of the coefficient of longi tudinal dispersion during the initial period before cross-sectional mixing is complete. The principal conclusions to emerge from this study are: 1. Numerical dispersion can be eliminated from the finite difference solution of mass transport equations by solving the equations along the advective characteris tics. • 2. The stability requirements of explicit finite difference schemes have been shown to be related to the speed of information propagation. Advective instabili ties are eliminated by solving the mass transport equa tion along the advective characteristics. 3. Existing theories have been shown to apparent ly underestimate the significance of secondary currents on the lateral mixing process. Secondary currents have been explained in terms of the generation and advection of vor-ticity, and estimated values show good agreement with limited field data. Revised estimates of lateral disper sion indicate significantly faster lateral mixing. 4. If the mass transport equation is solved along the advective characteristics, the time dependent behaviour of the coefficient of longitudinal dispersion during the initial period before cross-sectional mixing is complete can be taken into account. Also, a one-dimensional equa tion can be used to estimate lateral concentration profiles. 5. The effect of the tide has been shown to intro duce "spikes" into the concentration profile along the estu ary. These spikes are caused by variation in the initial dilution of a discharged effluent and multiple dosing due to tidal flow reversal. 6. Significantly more work and resources are in volved in the development and application of a tidally varying mass transport model than for a tidally averaged model. The time required for development, the amount of field data required for verification and the amount of computer time required to analyze an identical situation is estimated to be an order of magnitude greater for the tidally varying model as compared to the tidally averaged model. 94 7. The tidally averaged concentrations were not found to be a good indication of the tidally vary ing response of the Fraser River Estuary. The peak cross-sectionally averaged concentration predicted by the tidally varying model was found to be from one to 10 times greater than the tidally averaged values. The conclusions are now discussed in detail. In this study two mass transport models and a hydrodynamic model were developed and applied to the Fraser River Estuary, the hydrodynamic model being used to predict the temporal variation in the parameters of the tidally varying mass transport model. Finite difference methods have been used to solve the -underlying partial differential equations of all three models. Impor tant aspects of this type of solution are the problems of numerical dis persion and stability. There has been some confusion in the literature over the origin and means of controlling numerical dispersion. In this study, numerical dispersion has been shown to result from solving the mass transport equation over a fixed space grid rather than along the more fund amental advective characteristics; if the equation is solved along the ad vective characteristics, numerical dispersion is eliminated. As applied, Thomann's solution is steady state and has no numerical dispersion because it is independent of time. In addition to the problem of numerical disper sion, there is also the problem of stability in explicit finite difference schemes. In Appendix C this is seen to be related to the speed of infor mation propagation along the characteristics of the respective partial dif ferential equations. Stability requirements govern* the relative size of the space and time increments in the hydrodynamic equations and in the dispersion step of the tidally varying mass transport equation. There is no stability requirement for Thomann's solution as it is independent of time. However, 95 there is a limit on the relative magnitudes of the tidally averaged advec tive and dispersive transport processes, and if violated, the concentration in the offending segment will become negative. There are a number of advantages to solving the tidally varying mass transport equation along the advective characteristics; it results in a direct simulation of the advective transport occurring in the estuary, a useful separation of the advective and dispersive transport processes is achieved, and the position of each separate "effluent parcel" and the time that it has spent in the estuary is known. This last piece of infor mation has been used to account for an assumed > time-dependent increase in the coefficient of longitudinal dispersion in the initial period before cross-sectional mixing is complete, and also allows the tidally varying mass transport model to be adapted to predict lateral concentration profiles. A most useful feature of solving the equation along thej-advective charac teristics is that additional moving points can be placed on the character istics to more accurately define regions of rapid variation in concentra tion, such as occur at effluent outfalls at times of slackwater, and un necessary moving points can be removed from regions of slow variation. The concentration spikes generated by the effects of variation in initial dilu tion and multiple dosing are initially very sharp and only some 500 - 800 feet wide at the base. To achieve adequate resolution of such a spike with a fixed grid solution would require a space grid approximately 10 times finer than the standard 5,000 foot space grid used in this study. The advantages of solving the mass transport equation along the advective characteristics must be balanced against the somewhat untidy "book-keeping"of solution .results inherent to characteristics methods of solution. 95a In the tidally varying model of the Fraser River Estuary, this book-keeping is complicated by the muIti-channeled nature of the estuary and the inter actions of these channels at their junctions. But while awkward, the book keeping of results was not excessively difficult. Of the three models used in this study, only the hydrodynamic model has been verified with any degree of thoroughness. Because the tidally averaged dispersion coefficients have no real physical meaning, a tidally averaged model is forced to reproduce measured field results rather than being rigorously verified. Lack of field data prevented this, and also pre vented the complete verification of the tidally varying mass transport model. However, all three models are sufficiently flexible to simulate the range of flow and tidal conditions of the Fraser River Estuary and can be adjusted to fit field results when available. Essentially, the effect of the tides is to cause spikes in the concentration profile along the estuary. In a one-dimensional model, the tidally varying flows cause a variation in the initial dilution of a dis charged effluent, the concentration being greatest at times of slackwater. This generates spikes in the concentration profile along the estuary. The effects of tidal flow reversal result in certain slugs of water being dosed with effluent several times, and this also generates concentration spikes. Because of the asymmetric nature of the tides there is a period of weak flood and ebb flows once in each double tidal cycle. During this time there are three slackwaters and the velocities are low, and the effects of variation in initial dilution and multiple dosing interact to form a compound spike 96 whose concentration is significantly greater than that of the component spikes. After a spike has been generated, its concentration is reduced by the effects of lateral and longitudinal dispersion. The initial longitudinal dispersion after a "parcel of effluent" has been discharged into the estuary, is due principally to the effects of vertical velocity gradients. However, when the effluent is mixed over the cross-section, lateral velocity gradients dominate its longitudinal dispersion. To .account for this effect, it was assumed that the coefficient of longitudinal dis persion increased between these two extremes as described in Appendix F. Also, the coefficient of longitudinal dispersion was assumed to vary directly as the absolute velocity. While simplistic, these variations are thought to be a reasonable approximation of what occurs in the estuary. The peak tidally varying concentration at the point of effluent discharge was estimated to be from one to 10 times: greater than the value predicted from the tidally averaged model. Because the effluent is not uniformly distributed over the cross-section, the peak lateral concentration will be much greater than these values. After a spike has been in the estuary for several tidal cycles, its concentration has been greatly reduced by the longitudinal dispersion process. In the estuary with relatively short residence times such as the Fraser, lateral mixing will be of considerable importance. For example, in the lower reaches of the estuary, the residence time is only 2-4 tide cycles, and effluent may not be completely mixed over the cross-section even when it leaves the estuary. In the initial period before cross-sectional mixing is complete,the peak lateral concentration is significantly greater 97 than the cross-sectional averaged value predicted by the tidally varying model. However, the model can be adapted to predict lateral concentration profiles as was discussed previously. If cross-sectional mixing is not com plete at a diverging junction, the effluent may not be advected through the junction according to the simple flow balance of the mass transport models. To account for this effect in a model would require considerable field data. Because of the influence of secondary currents, cross-sectional mixing is thought to be relatively rapid in the Fraser River Estuary. Secondary flows have been tentatively explained in terms of the generation and advection of vorticity, and on the basis of limited float studies, show good agreement with values measured in the estuary. To sum up, the movement of water through the Fraser River Estuary is a complex phenomenon that is affected by the tides, the freshwater dis charge, the various channels and junctions of the estuary and the presence of Pitt Lake. The advective and dispersive transport processes that distribute an effluent throughout the- water mass of the estuary are similarly complex, and in particular are significantly affected by the tides. A tidally aver aged model, while simpler to develop and apply to the estuary, does not give a good indication of the tidally varying concentrations. This study does not pretend to solve the total problem of calculating effluent concentrations in an estuary as complex as the Fraser. Rather, it has concentrated on deve loping stable mathematical models free from numerical dispersion to allow a preliminary investigation of the significance of tidal effects. In addi tion an assessment is made of the importance of lateral dispersion, and it is seen that with simple modifications, the tidally varying model can at least partially account for the effects of lateral dispersion. REFERENCES Baines, W. D. [1952]. Water Surface Elevations and Tidal Discharges in the Fraser River Estuary, January 23 and 24, 1952, Report No. MH-32, National Research Council of Canada. Baines, w. D. [1953]. Survey of Tidal Effects on Flow in the Fraser River Estuary, June 10 and 11, 1952, Report No. MH-40, National Research Council of Canada. Bella, D. A. [1968]. "Solution of Estuary Problems and Network Programs," Discussion, Proc. A.S.C.E., J. Sanit. Eng. Div., 94 (SA1):180-181. Bella, D. A. and Greriney [1970]. "Finite-Difference Convection Errors," Proc. A.S.C.E., J. Sanit. Eng. Div., 96(SA6): 1361-1375. Bella, D. A. and Dobbins, W. E. [1968]. "Difference Modelling of Stream Pollution," Proc. A.S.C.E., J. Sanit. Eng. Div., 94'(SA5): 995-1016. Benedict, A. H., Hall, K. J. and Koch, F. A. [1973]. A Preliminary Water Quality Survey of the Lower Fraser River System, Technical Report No. 2, Westwater Research Center, University of British Columbia. Bennet, J. P. [1971]. "Convolution Approach to the Solution for the Dissolved Oxygen Balance in a Stream," Water Resources Res., 7(3): 580-590. Bird, R. B., Stewart, w. E., and Lightfoot, E. N. [I960]. Transport Pheno mena, (New York: John Wiley & Sons, 1960). Callaway, R. J. [1971]. "Application of Some Numerical Models to Pacific Northwest Estuaries," Proceedings, 1971 Technical Conference on Estuaries of the Pacific Northwest, Circular No. 42, Engineering Experiment Station, Oregon State University, Corvallis. Chow, V. T. [1959]. Open Channel Hydraulics, (New York: McGraw-Hill, 1959). Chow, v. T. [1964]. Handbook of Applied Hydrology, (Ed.) V. T. chow. (New York: McGraw-Hill, 1964), Section 8. 98 99 Courant, R. and Hilbert, D. [1962].', Methods of Mathematical Physics (2 vols; New York: Interscience, 1962), Vol. II. Custer, S. W. and Krutchkoff, R. G. [1969] . "Stochastic Model for BOD and DO in Estuaries," Proc. A.S.C.E., J. Sanit. Eng. Div., 95(SA5): 865-885. Diachishin, A. N. [1963], "Waste Disposal in Tidal Waters," Proc. A.S.C.E., J. Sanit. Eng. Div.3 89(SA4): 23-43. Di.Torb, D. M. and O'Connor, D. J. [1968]. "The Distribution of Dissolved Oxygen in a Stream with Time Varying Velocity," Water Resources Res., 4(3): 639-646. Di Toro, D. M. [1969]. "Stream Equations and The Method of Characteristics," Proc. A.S.C.E., J. Sanit. Eng. Div., 95(SA4): 699-703. Dornhelm, R. B. and Wodlhiser, D. A. [1968]. "Digital Simulation of Estuarine Water Quality," Water Resources Res., 4(6): 1317-1328. Dronkers, J. J. [1964]. Tidal Computations in Rivers and Coastal Waters, (Amsterdam: North-Holland Publishing Co., 1964). Dronkers, J. J. [1969]. "Tidal Computations for Rivers, Coastal Areas and Seas," Proc. A.S.C.E., J. Hyd. Div., 95(HY1): 29-77. Einstein, H. A. [1972]. Sedimentation, edited and published by Hsieh Wen Shen, Department of Civil Engineering, Colorado State University, Colorado. Elder, J. W. [1959]. "The Dispersion of Marked Fluid in Turbulent Shear Flow," J. Fluid Mech., 5: 544-560. Fischer, H. B. [1966a]. "A Note on the One-Dimensional Dispersion Model," Air Wat. Poll. Int. J., 10(6/7): 443-452. Fischer, H. B. [1966b]. Longitudinal Dispersion in Laboratory and Natural Streams, Report No. KH-R-12, Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California. 100 Fischer, H. B. [1967]. "The Mechanics of Dispersion in Natural Streams," Proceedings: ASCE Environmental Engineering Conference, Dallas, Texas, February 6-9, 1967. Fischer, H. B. [1969a], "The Effects of Bends on Dispersion in Streams," Water Resources Res.t 5(2): 496-506. Fischer, H. B. [1969b]. "Cross-Sectional Time Scales and Dispersion in Estuaries," Paper C. 19, Proceedings Volume 3 (Subject C), International Association for Hydraulic Research,, Thirteenth Congress, Kyoto, Japan, 31 August - 5 September, 1969. Fischer, H. B. [1970]. AlMethod for Predicting Pollutant Transport in Tidal Waters, Water Resources Center Contribution No. 132, Hydraulic Engineering Laboratory, College of Engineering, University of California, Berkeley. Fischer, H. B. and Holley, E. R. [1971], "Analysis of the Use of Distorted Hydraulic Models for Dispersion Studies," Water Resources Res., 7(1): 46-51. Fox, D. G. [1971]. "Finite Difference Convection Errors," Discussion, Proc. A.S.C.E., J. Sanit. Eng. Div., 97(SA5): 766-771. Gardiner, A. 0., Peaceman, D. W. and Pozzi, Jr., A. L. [1964]. "Numerical Calculation of Multi-Dimensional Miscible Displacement by the Method of Characteristics," J. Soc. Pet. Eng., 4(1): 26-36. Goldie, C. A. [1967]. Pollution and the Fraser: Report I Dealing with Pre liminary Investigations of Waste Disposal to the Lower Fraser River," prepared for the Pollution Control Board, Victoria, British Columbia. Guymon, G. L. [1970], "A Finite Element Solution of the One-Dimensional Diffusion-Advection Equation," Water Resources Res., 6(1): 204-210. Harder, J. A. [1971]. "Solution Techniques: 1; Analogue and Hybrid Tech niques," Chapter 6 of Estuarine Modelling: An Assessment, Ed. G. H. Ward and W. H. Espey, Jr. (listed under editors). 1G1 Harleman, D. R. F. [1965]. "The Significance of Longitudinal Dispersion in the Analysis of Pollution in Estuaries," Advances in Water Pollution Research, Ed. 0. Jaag (3 vols; Oxford: Pergamon Press, 1965), I, pp. 279-290. Harleman, D. R. F., Lee, C. H. and Hall, L. C. [1968], "Numerical Studies of Unsteady Dispersion in Estuaries," Proc. A.S.C.E., J. Sanit. Eng. Div., 94(SA5): 897-911. Harris, E. K. [1963]. "A New Statistical Approach to the One-Dimensional Diffusion Model," Air Wat. Poll. Int. J.,7(8): 799-812. Harris, E. K. [1964]. "Note on the One-Dimensional Diffusion Model and the Log-Normal Distribution," Air Wat. Poll. Int. J., 8(5): 421-423. Henderson, F. M. [1966]. Open Channel Flow (New York: The MacMillan Co., 1966). Hinze, J. O. [1959]. Turbulence: An Introduction to its Mechanism and Theory (New York: McGraw-Hill, 1959). Holley, E. R. [1969a]. "Unified View of Diffusion and Dispersion," Proc. A.S.C.E., J. Eyd. Div., 95(HY2): 621-631. Holley, E. R. [1969b]. "Difference Modelling and Stream Pollution," Proc. A.S.C.E., J. Sanit. Eng. Div., 95(SA5): 968-972. Holley, E. R., Harleman, D. R. F. and Fischer, H. B. [1970]. "Dispersion in Homogenous Estuary Flow," Proc. A.S.C.E., J. Hydv Div., 96(HY8).: 1691-1709. Kent, R. [1960]. "Diffusion in a Sectionally Homogenous Estuary," Proc. A.S.C.E., J. Sanit. Eng. Div., 86(SA2): 15-47. Koivo, A. J. and Phillips, G. R. [1971]. "Identification of Mathematical Models for DO and BOD Concentrations in Polluted Streams From Noise Corrupted Measurements," Water Resources Res., 7(4): 853-862. 102 Lager, J. A. and Tchobanoglous, G. [1968]. "Effluent Disposal in South San Francisco Bay," Proc. A.S.C.E., J. Sanit. Eng. Div., 94(SA2): 213-236. Leeds, J. V. [1967]. "Accuracy of Discrete Models Used to Predict Estuary Pollution," Water Resources Res., 3(2): 481-490. Leads, J. V. and Bybee, H. H... [1967]. "Solution of Estuary Problems and Network Programs," Proc. A.S.C.E., J. Sanit. Eng. Div., 93*'(SA3) : 29-36. Leendertse, J. J. [1971a]. "Environmental Simulation as a Tool in a Marine Waste Disposal Study of Jamaica Bay," Advances in Water Pollution Research, Ed. S. H. Jenkins (2 vols; Great Britain: Pergamon Press, 1971), U, Sec tion III, pp. 4/1-4/7. Leendertse, J. J. [1971b], "Solution Techniques: 2; Finite Differences," Chapter 6 of Estuarine Modelling: An Assessment, Ed. G. H. Ward and W. H. Espey, Jr. (listed under editors). Li, W. H. [1962]. "Unsteady Oxygen Sag in a Stream," Proc. A.S.C.E., J. Sanit. Eng. Div., 88(SA3): 75-85. Li, W. H. [1972]. "Effects of Dispersion on DO Sag in Uniform Flow," Proc. A.s.C.E., J. Sanit. Eng. Div., 98(SA1): 169-182. Loucks, D. P. [1967]. "Risk Evaluation in Sewage Treatment Plant Design," Proc. A.S.C.E., J. Sanit. Eng. Div., 93(SA1): 25-39. Loucks, D. P. and Lynn, W. R. [1966]. "Probablistic Models for Predicting Quality," Water Resources Res., 2(3): 593-605. NalTuswami,; M.^, .Longenbaugh,; R. A. and Sunada, D." K. [1972]. "Finite Element Method for the Hydrodynamic Dispersion Equation with Mixed Partial Deriva tives," Water Resources Res., 8(5): 1247-1250. O'Connell, R. L. and Walter, C. M. [1963]. "Hydraulic Model Tests of Estuarial Waste Dispersion," Proc. A.S.C.E., J. Sanit. Eng. Div., 89(SA1): 51-65. O'Connor, D. J. [I960]. "Oxygen Balance of an Estuary," Proc, A.S.C.E., J. Sanit. Eng. Div., 86(SA3): 35-55. 103 O'Connor, D. J. [1962]. "Organic Pollution of New York Harbour - Theoreti cal Considerations," J. Water Poll. Control Fed.3 34(9): 905-919. O'Connor, D. J. [1965]. "Estuarine Distribution of Non-Conservative Sub stances," Proc. A.S.C.E., J. Sanit. Eng. Div., 91(SA1): 23-42. O'Connor, D. J. [1966]. "An Analysis of the Dissolved Oxygen Distribution in the East River," J. Water Poll. Control Fed., 38(11): 1813-1830. O'Connor, D. J. [1967]. "The Temporal and Spatial Distributions of Dissolved Oxygen in Streams," Water Resources Res.3 3(1): 65-79. O'Connor, D. J., St. John, J.P. and Di Toro, D. M. [1968]. "Water Quality Analysis of the Delaware Estuary," Proc. A.S.C.E., J. Sanit. Eng. Div.3 94(SA6): 1225-1252. O'Connor, D. J. and Di Toro, D. M. [1970]. "Photosynthesis and Oxygen Balance in Streams," Proc. A.S.C.E., J. Sanit. Eng. Div.3 96(SA2): 547-571. Odd, N. V. M. [1971]. "Practical Application of Mathematical Model Techniques to Predict Effects of Engineering Works on Tidal Flow in Estuaries With A Large Tidal Range," Proc. Inst. Civ. Engrs.3 50(December): 507-517. Oster, C. A., Sannichsen, J. C. and Jaske, R. T. [1970]. "Numerical Solution to the Convective-Diffusion Equation," Water Resources Res.3 6(6): 1746-1752. Pence, G. D. Jr., Jeglic, G. D. and Thomann, R. V. [1968], "Time Varying Dissolved Oxygen Model," Proc. A.S.C.E., J. Sanit. Eng. Div.3 94(SA2): 381-402. Pinder, G. F. and Cooper, Jr. H. H. [1970]. "A Numerical Technique for Calculating the Transient Position of the Saltwater Front," Water Resources Res., 6(3): 875-882. Pinder, G. F. and Friend, G. 0. [1972]. "Application of Galerkin's Procedure to Aquifier Analysis," Water Resources Res., 8(1): 108-120. Preddy, W. S. and Webber, B. [1963]. "The Calculation of Pollution of the Thames Estuary by a Theory of Quantized Mixing," Air Wat. Poll. Int. J., 7(8): 829-843. 104 Pretious, E. S. [1972]. Downstream Sedimentation Effects of Dams on the Fraser River, B.C., Water Resource Series Report No. 6, Department of Civil Engineering, University of British Columbia. Price, H. S., Cavendish, J. C. and Varga, R. S. [1968], "Numerical Methods of Higher-Order Accuracy for Diffusion-Convection Equations," J. Soc. Pet. Eng., 8(3): 293-303. Pritchard, D. W. [1969]. "Dispersion and Flushing of Pollutants in Estuaries," Proc. A.S.C.E., J. Hyd. Div., 95(HY1): 115-124. Prych, E. A. and Chidley, T. R. E. [1969]. "Numerical Studies of Unsteady Dispersion in Estuaries," Discussion, Proc. A.S.C.E., J. Sanit. Eng. Div., 95(SA5): 959-964. Quick, M. C. [1973], "The Mechanics of Streamflow Meandering," paper submitted for publication. Rennerfelt, J. G. W. [1963]. "An Analogue Computer for Oxygen Sag Calcu lations," Air Wat. Poll. Int. J., 7(4/5): 367-372. Richtmyer, R. D. and Morton, K. w. [1967]. Difference Methods for Initial Value Problems. 2nd Ed. (New York: Interscience, 1967). Sawyer, W. W. [1966]. A Path to Modem Mathematics (London: Penguin Books Ltd., 1966). Schofield, W. R. and Krutchkoff, R. G. [1972]. "Stochastic Model for BOD and OD in Segmented Estuary," Proc. A.S.C.E., J. Sanit. Eng. Div., 98 (SA5): 745-762. Simmons, H. B. [1960]. "Application and Limitation of Estuary Models in Pollution Analyses," Proceedings: First International Conference on Waste Disposal in the Marine Environment, Ed. E. A. Pearson (Oxford: Pergamon Press, 1960), pp. 540-546. Streeter, H. w. and Phelps, E. B. [1925]. A Study of the Pollution and Natural Purification of the Ohio River, Public Health Bulletin 146, U.S. Public Health Service, Washington, D.C. 105 Taylor, G. I. [1954] . "The Dispersion-, of Matter in Turbulent Flow Through a Pipe," Proa. Roy. Soc. London, 223 (Series A): 446-468. Thayer, R. P. and Krutchkoff, R. G. [1967]. "Stochastic Model for BOD and DO in Streams," Proc. A.S.C.E., J. Sanit. Eng. Div., 93(SA3): 59-72. Thomann, R. V. [1963]. "Mathematical Model for Dissolved Oxygen," Proc. A.S.C.E., J. Sanit. Eng. Div., 89(SA5): 1-30. Thomann, R. V. [1965]. "Recent Results From A Mathematical Model of Water Pollution Control in the Delaware Estuary," Water Resources Res., 1(3);: 349-359. Thomann, R. V. [1971]. Systems Analysis and Water Quality Management (New York: Environmental Science Services Division, Environmental Research and Applications, 1971). Waldichuk, M. Market, J. R., and Meikle, J. H. [1968]. Fraser River Estuary, Burrard Inlet, HoweSound and Malaspina Strait Physical and Chemical Ocean ographic Data, 1957-1966: Volume 1, September 1957 to February 1962, Manuscript Report Series No. 939, Fisheries Research Board of Canada. Ward, G. H. and Espey, w. H. [1971]. Estuarine Modelling: An Assessment, Ed. G. H. Ward and W. H. Espey, Report No. 16070 DZV 02/71, Water Pollution Research Series, Environmental Protection Agency, U.S.A. Ward, P. R. B. [1972]. Transverse Dispersion of Pollutants in Oscillatory Open-Channel Flow, Report No. WHM-2, Department of Civil Engineering, University of California, Berkeley. Ward, P. R. B. [1973]. "Prediction of Mixing Lengths for River Flow Gauging," Proc. A.S.C.E., J. Hyd. Div., 99(HY7): 1069-1081. Worley, J. L., Burgess, F. J. and Towne, W. W. [1965]. "Identification of Low-Flow Requirements for Water Quality Control by Computer Techniques," J. Water Poll. Control Fed., 37(5): 659-673. APPENDICES APPENDIX A DERIVATION OF THE ONE-DIMENSIONAL MASS TRANSPORT EQUATION A. 1 GENERAL The mass transport processes of advection, turbulent diffusion and molecular diffusion transport and distribute any dissolved substance throughout the water mass of a river or estuary. (Advection is also re ferred to as convection or forced convection in the literature). At the time and space scales of interest, the contribution to mass transport by molecular diffusion is insignificant compared to the other processes, and is not considered further. The inability of the advective or bulk velocity to describe the distribution of velocity over a cross-section of the river or estuary makes, it necessary to postulate another transport mechanism, namely dispersion. The mechanisms of advective and dispersive transport are described in Sections A.2 and A.3. The one-dimensional mass transport equation is derived by taking a mass balance of dissolved substance over an element cross-sectional slice of the estuary. Dissolved substance is transported into and out of this volume by the processes of advection and dispersion, and these processes, together with any source-sink effects that the substance undergoes within the water mass, determine its concentration within the elemental volume. The only spatial variation that is admitted in the one-dimensional mass transport equation is a longitudinal variation along the estuary. The concentration, longitudinal velocity, depth and dispersion coefficient are 106 107 all assigned their average cross-sectional values (c, u, y and E respec tively) . For the sake of brevity, the average cross-sectional value of any parameter will be referred to as its "mean" value. The high fre quency turbulent fluctuations are assumed to have been averaged out of the mean values u and c. Figure A.l shows an elemental cross-sectional slice of an estuary. The direction of the co-ordinate axes x, y and z are as shown, the positive x axis pointing seawards and the positive y axis vertically down. Note that the co-ordinate axes are Eulerian, or fixed in space. The flow is assumed to be both unsteady and non-uniform so that c, u, y, E and the cross-sectional area A vary both with time (t) and longitudinal distance (x). The width of estuary b is assumed to vary with x only. It is assumed that the longitudinal variation in all these parameters is con tinuous. A.2 LONGITUDINAL ADVECTIVE TRANSPORT The distribution of longitudinal velocity in an estuary is three-dimensional. Variations in the longitudinal direction occur due to changes in the cross-sectional area of the estuary, and over any cross-section, variations in the lateral and vertical directions occur due to the frictional influence of neighbouring fluid layers with each other and the solid bed. In a river or estuary, the lateral distribution of longi tudinal velocity is modified by the presence of bends, and in an estuary, the vertical distribution of longitudinal velocity is modified by the additional influence of salinity. 108 Figure A.l Elemental Cross-Sectional Slice of A River or Estuary 109 In one-dimensional mass transport models, it is assumed that a cross-sectional slice or slug of water moves along the river or estuary at the mean cross-sectional velocity u. The transport of dissolved sub stance in this slug of water is termed longitudinal advection. The advective transport in a river is always in the downstream direction. In an estuary, the net advective transport over a tidal cycle is in the downstream direction because of the effect of freshwater inflow. However, because of flow reversal due to tidal effects, the advective transport will be in the upstream direction over portion of the tidal cycle. The mass of dissolved substance (M) advected through any cross-section of a river or estuary in time 6t is given by M = uAc 6t, and thus the net mass (M \ remaining in the elemental volume of Figure net A.l due to advective transport during time 6t is g - -M = - r—[uAc}6x6t (A.l) net 9x A.3 LONGITUDINAL DISPERSIVE TRANSPORT Advection accounts for the transport of dissolved substance in the longitudinal direction due to the assumed uniform distribution of the velocity over the cross-section, u.,,The vertical and lateral velocity gradients that exist over any cross-section result in small "parcels" of water preceding and lagging the line of advective advance, as defined by u. This is illustrated for a vertical velocity profile typical of a river 110 Vertical Distribution of Velocity UAt V7 A o —~ o o -^Line source Y position after o Time At 0 >v 0 > a 0 o i 0 0 ° o ° 1 ° o o o o o o o o o o A Effect of Diffusion f^^^--Advection Due to Velocity Gradient Figure A.2 Dispersive Effects of Vertical Velocity Gradients Ill in Figure A.2. It is seen that a slug input from a line source is not advected downstream as a slug load, but is spread around the line of mean advance. The effect of this is to reduce the peak concentration and flatten the concentration gradients, as predicted on the basis of advection alone. Superimposed on the effects of velocity gradients is the process of turbulent diffusion. The mass transport associated with the turbulent velocity fluctuations results in a further flattening of the concentration gradients. The effect df turbulent diffusion over the vertical section "AA" is also illustrated in Figure A.2. The combined effects of velocity gradients and turbulent diffusion in spreading the dissolved substance around the line of advective advance (u) is termed longitudinal dispersion. Unless . otherwise qualified, the terms "advection" and "dispersion" will be taken to mean ad vection and dispersion in the longitudinal direction. It is necessary to represent the dispersion process mathematically if it is to be included in the mass transport equation. Holley [19 69a]dis cussed the underlying similarity between the mechanisms of molecular diffu sion, turbulent diffusion and dispersion. Both molecular diffusion [Bird, et al.s 1960] and turbulent diffusion in homogeneous isotropic free turbu lence [Hinze, 1959] can be represented mathematically by a Fickian diffusion equation. Taylor [1954] showed that after a suitable time had elapsed, the longitudinal dispersion in flow through a pipe can also be represented by a one-dimensional Fickian diffusion equation. Elder [1959] found a similar result to hold in model experiments of open channel flow. Fischer [1966a] illustrated that in the initial non-Fickian period following the release of tracer, the effects of velocity gradients outweighed the effects of turbu lent diffusion, and led to a skewed distribution of mean concentration in the longitudinal direction, as is illustrated in Figure A.2. According to Fick's Law, the distribution of concentration should be Gaussian, but Fischer [1966a] has shown that this only occurs after cross-sectional mixing has reduced the spatial variation of concentration over the cross-section to a value much smaller than the mean cross-sectional value. In deriving and applying the one-dimensional mass transport equation, the dispersion process in its entirety is assumed to be Fickian. In actual fact, the longitudinal dispersion of effluent in a river or estuary is a complex phenomenom that is initially controlled by the effects of vertical velocity gradients. As turbulent diffusion and secondary flows distribute the effluent mass over the cross-section, the lateral velocity gradients exert an increasing influence on the longitudi nal dispersion process, The lateral mixing due to turbulent diffusion and secondary flows is discussed in Appendix E and an estimate of the coeffi cients of longitudinal dispersion of the estuary is given in Appendix F. (The effects of vertical and lateral velocity gradients are also discussed in Appendix F). When the cross-sectional mixing is essentially complete (and the dispersion process is Fickian) the effects of lateral velocity gradients dominate the longitudinal dispersion and are of the order of 20 - 100 times greater than the effects of vertical velocity gradients [Fischer, 1966a]. In an estuary, the dispersion process is modified by the additional influence of salinity and tidal effects. The influence of tidal oscillations on the dispersion coefficient has been investigated by Holley et al. [1970] , and their method has been used in estimating the dispersion coefficients for the Fraser River Estuary in Appendix F. According to Fick's law of diffusion, the net mass of dissolved substance transported through unit cross-sectional area in time 6t is given by [Bird et al., 1960] M = -E—fit (A. 2) 3x 113 the minus sign arising because the net transport is in the direction of decreasing concentration. E is the dispersion coefficient, and as is expected from the previous discussions of the dispersion process, depends on the distribution of velocity over the cross-section and the coefficient of turbulent diffusion (see, for example, Taylor [1954]; Fischer [1966a], [1967] and Holley et al. [1970]. According to equation (A.2), the net mass remaining in the elemental volume due to dispersive transport during time 6t is M . = ^-tAE^KxSt. (A.3) net 3x 3x A.4 SOURCE-SINK EFFECTS The source-sink terms account for any processes that result in the production or removal of the dissolved substance from the water body. In the BOD-DO system, source processes include reoxygenation by atmospheric oxygen dissolving in the water and the oxygen produced photosynthetically by aquatic plants; the sink processes include oxygen uptake by bacteria breaking; down the substances exerting the BOD and oxygen uptake for aquatic plant respiration. The source-sink processes that a substance undergoes in the aqua-•tic, environment are specific to the particular substance, and depend on the physical, chemical and biological processes that the substance undergoes in the particular aquatic environment. Consequently, source-sink effects are included in the mass transport equation dh a general manner. Suppose that the particular substance is undergoing n processes that result in its produc tion or removal from the water body. If the rate of production per unit 114 volume of water due to the i process is given by S^, removal processes being regarded as negative production, then the net, mass of substance pro duced in the elemental volume in time 6t is n M = A&x&tl S. (A.4) net . , 1 i=l A.5 THE ONE-DIMENSIONAL MASS TRANSPORT EQUATION The mass of dissolved substance contained in the elemental volume is M = A<5xc As a result of the advection, dispersion and source-sink effects that occur in time 6t, the change in mass of dissolved substance in the ele mental volume is 3 — 6M = T--{Ac}6x6t (A.S) dt Summing equations (A.l), (A.3) and (A.4) and substituting into equation (A.5) gives |-CAc}> 4-iAuc"} + |-[AE|£-} + A E S. (A.6) 9t 3x 3x 3x . , l i=l Expanding the first two derivatives and using the equation of conservation of fluid mass in the form ' 3 <• -i 3A N (Au) + — = 0, 3x dt equation (A.6) reduces to 3c 3c 1 3 f 3c •, . sr = -ur— +7-1 <AE;r— } + I S., (A.7) 3t 3x A 3x 3x . , l 1=1 115 the one-dimensional mass transport equation for unsteady non-uniform flow. The terms on the right-hand side of Equation (A.7) shall be re ferred to respectively as the advective, the dispersive and the source-sink terms. A.6 THE LAGRANGIAN FORM OF THE ONE-DIMENSIONAL MASS TRANSPORT EQUATION Equation (A.7) was derived relative to an Eulerian co-ordinate system (xyz) that is fixed in space. To transform the equations relative to a co-ordinate system that is advected along at the mean velocity u (a Langrangian co-ordinate system), it is necessary to transform the Eulerian co-ordinates according to §- = u (A.S) dt Relative to the curve (A.8), the time rate of change of concentration is given by (see Section Cl) dc _ 3c_ -3c dt " 3t + U3x and thus, relative to Langrangian co-ordinates, Equation (A.7) becomes £ = | T-CAE& + Z S. (A.9) dt A 3x 3x ,i x=l Equation (A.8) is the characteristic curve of advective information propagation, and is discussed further with regard to numerical dispersion and stability in Appendix C. Equation (A.7) and Equations (A.8) and (A.9) repre sent the same situation but from different viewpoints; in Equation (A.7) the "observer" is stationed on the bank of the river or estuary and observes changes 116 in the concentration as the water flows past him, whereas in the latter two equations, the observer is advected along with the water mass [Equation (A.8)] and observes changes in the water mass around him [Equation (A.9)]. If dispersive transport is ignored, Equation (A.8) reduces to dt i-1 1 which is the original one-dimensional mass transport equation developed by Streeter and Phelps [1925] to predict the effect of substances exerting a biological oxygen demand (BOD) on the concentration of dissolved oxygen (DO) in a river. In their model, n was equal to two; the only source of DO was reaeration through the water surface and the only sink was the satisfaction of the BOD. APPENDIX B A DESCRIPTION OF THE FRASER RIVER ESTUARY B.l GENERAL The Fraser River is the largest and most important river in British Columbia. Its drainage area of 90,000 square miles is shown in Figure B.l and occupies almost one-quarter of the Province. The river rises in. the Rocky Mountains near Jasper and flows north-west for some 330 miles. • North of Prince George the river changes direction and flows south some 400 miles to Hope. At Hope the river emerges from the Coast Mountains and flows westward through the Fraser Valley for almost 100 miles to enter the Strait of Georgia at Vancouver. It is this latter section of the river, from Hope to Vancouver, that is of direct concern to this thesis, and dis cussion is now restricted to it, B.2 THE LOWER FRASER RIVER SYSTEM The Lower Fraser River from Hope to Vancouver is shown in Figure B.2. A convenient way of describing this section of the river is in terms of the factors that affect the flow. These factors include: (1) the var ious channels of the river system; (2) the tides in the Strait of Georgia; (3) the river discharge; and (4) saltwater intrusion. B.2.1 Channels of the River System. The Fraser River deposits 20 - 25 million tons of sediment into the Strait of Georgia annually [Pretious, 1972], and has built a large delta with a seaward perimeter of some 30 miles. 117 00 The Fraser River From Hope to Vancouver 120 The delta is shown in Figure B.3 and it is seen that the river has cut a net work of channels through it. At New Westminster the river trifurcates into the Main Arm, which enters the Strait off Steveston; Annacis Channel, which flows around the back of Annacis Island to reconnect with the Main Arm; and the North Arm, which enters the Strait near Point Grey. Annacis Channel is not treated as a separate channel, but as part of the Main Arm. Around Sea Island, the North Arm divides to form the Middle Arm. North-west of Ladner, a complex of shallow channels and sloughs (backwater areas) interconnect with each other, the three most important channels being Ladner Reach, Sea Reach and Canoe Pass. This area?; drains through Canoe Pass to provide the Main Arm with another connection to the Strait of Georgia. The percentage of river discharge that enters the Strait through the Main Arm, the North Arm, the Middle Arm and Canoe Pass is approximately 85 per cent, 5 per cent, five per cent and five per cent respectively [Goldie, 1967-]. This is only a general indication of the percentage distribution of flows, the summer dis tribution of flows probably differing from the winter distribution. The section of river between New Westminster and Hope is referred to as the Main Stem. Of the various islands present in the Main Stem (see Figure B.2), only Douglas and Barns ton Islands form significant additional channels. The channels at the "back" of the other islands are treated as part of the Main Stem. Of the various major lakes that drain into the Main Stem (see Fig ure B.2), Pitt Lake is of particular significance. The surface elevation of the lake is affected both by the tides in the Strait of Georgia and the fresh water discharge of the Fraser (the discharge of Pitt River into the north end 121 Point Ladner Reach LADNER V;'v'::"';--':'"•-"•^•'.f' Canoe 'Pass' SCALE: 12 3 4 5 MILES LEGEND : ^Y::i:;; - Delta CANADA U , S . A Figure B.3 The Fraser River Delta 122 of the lake is negligible compared to the Fraser discharge. At low fresh water flows in the Fraser, the maximum tidal range of the lake is several feet, and as the lake has a surface area of some 25 square miles, the volume of water stored in the lake on the flood tide and released to the Fraser on the ebb tide is very large (see Table B.l). The large reverse delta at the entrance to Pitt Lake (see Figure B.2) is evidence of the large flows that enter the lake. The Pitt River - Pitt Lake system obviously consti tute an important channel of the lower Fraser,[/system and will be referred to as the Pitt System. The other lakes draining into the Main Stem are of sufficient elevation to be independent of tidal influence and freshwater discharge in the Fraser. Thus, the lower Fraser River system consists of the following seven principal channels: 1. Main Arm - Main Stem 2. North Arm 3. Middle Arm 4. Canoe Pass 5. Pitt System 6. Douglas Channel, and 7. Barnston Channel. To some extent the channel geometries are determined by the many .-•.. dykes, jetties and training walls that have been constructed in the lower Fraser system. Each Spring the Fraser Valley is subject to possible flooding during the freshet (time of high runoff due to snowmelt), and an extensive system of dykes has been constructed for flood protection. Ocean-going ships 123 use the Main Arm to enter and leave New Westminster, an overseas shipping terminal. To maintain the shipping channel, a number of dykes and training walls have been constructed to increase local scour. Cross-sectional areas, depths and widths were evaluated at the net work of stations shown in Figure B.4. The stations are 5,000 feet apart, except in the deeper water of Pitt Lake where the spacing has been increased to 15,000 feet. (This network of stations is used in the finite difference solutions of the hydrodynamic equation and the tidally varying and tidally averaged mass transport equations, as is discussed in Chapter 3). Soundings charts supplied by the Department of Public Works of Canada were used to deter mine the channel geometries. A typical cross-section of each of the seven channels is shown in Figure B.5. Widths and depths are recorded in feet, and depths are relative to local low water. Figures B.6 to B.8 show the var iation in depths, widths and cross-sectional areas along the five major channels. In an attempt to estimate the advective area (effective area of discharge), it was arbitrarily assumed that any depth of section less than 10 feet deep did not contribute to the advective area. This is illustrated for Station No. 35 in Figure B.5, where areas less than 10 feet deep are shown hatched. The storage width of section is AB, but the advective width is only CD. The advective area was used to determine the advective hydraulic radius (which is greater than the hydraulic radius determined from the gross area of section). From Figures B.6 to B.8, it is seen that the advective area is not much smaller than the gross area, but the advective width is often considerably smaller than the storage width. The river sections are essentially wide and shallow and the terms hydraulic radius and mean depth are now used synonymously. FIRST AND LAST STATIONS •• Main Arm - Main Stem I North Arm 101 Canoe Pass 11 9 Middle Arm 129 Pitt System 140 Douglas Channel 180 Barnston Channel 184 , North Arm 101 S J06 HQ, •Middle Arm ®jS "Main Arm *I28 Pitt Lake Junction Stations 5000* spacing 15,000' spacing •Canoe Pass Figure B.4 Network of Stations Used in the Numerical Solution of the Hydrodynamic and Mass Transport Equations 1000 2000 125 Stat. No. 35 ( Main Stem ) Stat. No. 109 ( North Arm) 1000 Stat. No. 134 ( Middle Arm ) 1000 Stat. No. 126 ( Canoe Pass) 1000 1— , Stat. No. 145 ( Pitt) Stat.No. l8l(Douglas Channel) Stat. No. I87( Barnston Channel) Figure B.5 Typical Channel Cross-Sections (Widths and Depths in Feet, Depths Relative to Local Low Water) 60,000 h 50,000 o 40,000 30,000 20,000 50 CD < CL OJ Q 10 3,000 sz 2,000 •D 1,000 0 Gross Values Advective Values DEPTHS AND AREAS Relative to local low water. 10 20 30 Stations 40 50 Figure B.6 Cross-Sectional Parameters of Main Arm - Main Stem 60 Figure B„7 Cross-Sectional Parameters of North Arm, Middle Arm and Canoe Pass DEPTHS AND AREAS Relative to local low water. Pitt River 40,000 i; 30,000 o 20,000 a> k. < 10,000 0 40 £ 30 20 I 0 0 3,000 - 2,000 CL a> Q •o ,000 o z o in o z Pitt Lake in m in o z o (0 o z co in ro o Z /A y. r^yjr§e_ Delta 1 \\ H- \ »\ -V — ' » / — / /^\\ / V » \ ~—v7\\ §W 'SN—\ % \ t \ % \ V* — i I / —H 1 1 — Gross Values Advective Values 2,000,000 ~ 1,500,000 ~ 1,000,000 < 500,000 0 2 10 ^ 200 ~ 190 180 CL Q 10,000 8,000 _ 6,000 -JZ 4,000 5 2,000 03 Figure B.8 Cross-Sectional Parameters of Pitt River and Pitt Lake 129 B.2.2 Tides in the Strait of Georgia. The Strait of Georgia is tidal, being connected to the Pacific Ocean by the Juan de Fuca Strait (see Figure B.l). The tides are of the mixed type characteristic of much of the coast of Northwest America. The tidal range at Steveston for mean and large tides is respectively 10 feet and 15 feet and typical tides at Steveston are shown in Figure B.9. The effects of the tides on the flows in the Fraser system depends both on the tidal range at Steveston and the river discharge. This is illus trated in Figure B.10 [after Baines, 1953] which shows the local low and high tide envelopes of the Main Arm and Main Stem for discharges of 27,000 and 250,000 cubic feet per second at Hope. The greater influence of the tide during low flow conditions is readily apparent, and it is seen that the upstream limit of tidal influence is around Chilliwack, some 60 miles from Steveston. The section of the Fraser system from the Strait of Georgia to Chilliwack will be referred to as the Fraser River Estuary, although it is noted that the Fraser is more properly designated as a tidal river [Callaway, 1971]. During low flow-high tide conditions, flow reversal occurs in the Fraser Estuary. The cubature study of discharges by Baines [1952] predicted flow reversal at Mission, some 50 miles upstream of Steveston. Downstream of Chilliwack there is a network of tide gauging stations, . some.equipped with continuous recorders and others only recording the maximum .and minimum levels during each 24 hour period. The position and type of gauge at each station is shown in Figure B.ll. • : B.2.3 River Discharges. The variation of t^g'meap^n^n.thly" "dis^* charge'~?at Hope, as calculated for the period 1912-1970;, is':sriowh;.-in Figure-''. B. 12. .'The .mean monthly flows vary from a summer maximum of-250,000""cfs during 0.0 3.0 6.0 9.0 12.0 15.0 18.0 21.0 24.0 Hours 6.0' 1 1 1 1 1 1 1 I 0.0 3.0 6.0 9.0 12.0 15.0 18.0 21.0 24.0 Hours Figure B.9 Typical Tides at Steveston Steveston New Mission Chilliwack Hope Westminster Figure B.J.0 Local Low and High Tide Envelopes [After Baines, 1953] Figure B.ll Tide Gauging Stations in the Fraser River Estuary to 133 280i 240 (/> "2001 w O O — I 601 o LL. 120 80 40 0' JAN. FEB. MAR. APR. MAY JUNE JULY AUG. SER OCT. NOV. DEC. Figure B.12 Mean Monthly Flows at Hope (1912-1970 Inclusive) 134 freshet 7:to a winter minimum of around 30,000 cfs. Recorded extremes at Hope are 536,000 cfs on May 31, 1948 and 12,000 cfs on January 8, 1916. The drainage area of the Fraser River below Hope is some 6,000 square miles, and between Hope and the Strait of Georgia various other rivers flow into the Fraser system, the most important being the Harrison River (see Figure B.2). The effect of this additional inflow is to increase the flow at New Westminster by some 15 per cent over the flows at Hope during the freshet, and up to 50 per cent during winter. It should be apparent that the Fraser River Estuary is somewhat unusual, being characterized by both high river discharges and large tidal effects. This is illustrated by a calculation of the tidal prism and total river discharge between the times of low-low-water and high-high-water for the tidal cycles of January 11 and June 16, 1964. These tidal cycles and the freshwater discharges are shown in Figure B.9. The volumes of the tidal prisms and river flows are shown in Table B".l. B.2.4 Saltwater Intrusion. The salinity of the Strait of Georgia is essentially some 30 parts per thousand (ppt), and under low river discharge conditions saltwater intrudes a considerable distance into the four channels that emerge from the delta. The saltwater intrusion into the Main Arm is illustrated in Figure B.13 [data from Waldichuk, et al. , 1968], which shows the salinity profile at four stations along the Main Arm. (Note: the station designation in Figure B.13 is that of Waldichuk, et al. 3 and should not be confused with the station designation in Figure B.4). The position of the stations, the tide at Steveston and the time of tide the observations were taken are also shown. The four observations closely occur around the same phase of tide and can be regarded as "simultaneous." The freshwater discharge TABLE B.l RIVER FLOW VOLUMES AND TIDAL PRISMS ON JANUARY 15 AND JUNE 16, 1964 (Volumes Calculated for Time Between Low-Low-Water and High-High-Water) DATE RIVER FLOW (cfs) TIDAL RANGE (feet) RIVER FLOW vol. (ft ) TIDAL PRISM (ft3) RIVER FLOW TIDAL PRISM JANUARY 15 53,500 10.4 12 x 108 66 x 108 0.18 JUNE 16 463,000 7.6 104 x 108 25 x 108 4.2 10 20 30 •° • -JL. 1 10 20 30 10 20 30 100 200 300L All depths in feet River Discharge - 57,600 cfs at Hope 10 20 30 —I 1 1 Stat. No. 3 Figure B.13 Salinity Profiles in the Main Arm on February 13-14, 1962 01 137 at Hope is 57,000 cfs and it is apparent that under these freshwater dis charge and tidal conditions, the Main Arm is highly stratified with the toe of the salt wedge somewhere between stations 2 and 3. Recent studies have involved monitoring the salinity profile con tinuously at three stations along the Main Arm. This was done for the two periods February 1-16 and March 16 - 30, 1973, during which the average discharge at Hope was 30,000 cfs ^and 34,000 cfs respectively. It was found that the effect of the tides was to move the wedge bodily up and down the Main Arm [D. O. Hodgins, private communication], the maximum excursion of the toe of the wedge probably being to somewhere around Annacis Island. Thus, for low river flows, the Main Arm of the Fraser Estuary is highly stratified with the salt wedge moving bodily up and down the estuary under the influence of the tides. During times of flow reversal, the wedge may move as far upstream as Annacis Island but during seaward discharge it is washed downstream past Steveston. For high river flows, the wedge prob ably does not penetrate'-past Steveston, if that far. Saltwater intrudes into the other channels of the delta, and similar situations probably occur there, although the saltwater movement may be modified by the slower velocities through these channels. Salinity profiles in these other channels show them to be typically stratified, but perhaps not quite as stratified as the Main Arm (probably due to the lower discharge velocities through these channels). APPENDIX C NUMERICAL DISPERSION AND STABILITY The dependent variable of a partial differential equation de fines a surface over the plane of the independent variables. In the plane of the independent variables there are various curves, or "characteristics" that describe the propagation of information through the system that the partial differential equation represents. ' The one-dimensional forms of the advective transport equation, the dispersive transport equation and the hydrodynamic equations are of interest to this thesis. The forms of their respective surfaces are briefly described and the equations of their respec tive characteristics are given. In fixed grid finite difference schemes, the problem of numerical dispersion arises from solving the one-dimensional advective transport equation (or the one-dimensional mass transport equation) relative to a fixed space grid, rather than along the advective characteris tics. In explicit finite difference schemes, the problem of stability also arises from solving the above three equations relative to a fixed grid rather than along their respective characteristics. Cl SURFACE GEOMETRY OF PARTIAL DIFFERENTIAL EQUATIONS The following partial differential equations are of interest 3c 3t (Cl) (C.2) 138 139 9u 9t *37 -?f= 9a u u c y (C.3) 9y_ 9t -9u Equation (Cl) describes one-dimensional advective transport and Equation (C.2) nent transport processes of the one-dimensional mass transport equation, and the sum of their effects is the one-dimensional mass transport equation with out the source-sink terms. The two coupled equations (C.3) are the hydrodynamic equations and are used to predict the temporal variations in the parameters u and A of the tidally varying mass transport equation. [Equations (C.3) are ob tained by substituting Equation (3.3) and A = by into Equations (3.1) and (3.2)]. Equations (Cl) and (C.2) define surfaces c(x,t) over the (x,t) plane and Equations (C.3) define a pair of coupled surfaces u(x,t) and y(x,t) over the (x,t) plane. Consider now the form of the surface defined by the advective transport Equation (Cl). Figure Cl shows the shape of this surface for a slug input into (a) steady uniform flow; (b) steady non-uniform flow along a river of contracting cross-section; and (c) the mixture of steady and oscillatory flow characteristic of an estuary (estuary flow). These sur faces are solutions to Eguation (Cl) for a slug input into the various types of flow. It is apparent from Figure Cl that there are various curves in the (x,t) plane that result in a considerable simplification of Equation (C.l) — the concentration is constant along each of the curves AB of Figure C.l while everywhere else in the plane it is zero. For any continuous, smoothly varying surface c(x,t), the time rate of change of c(x,t) along the vertical plane defined by the curve one-dimensional dispersive transport.'"These two equations represent the compo-140 c(x,t) Figure Cl Concentration Surfaces for the Advection of a Slug Load 141 'Sr = f(x,t) (G.4) at in the (x,t) plane is given by dc _ 3c 3c_ dx dt 3t 3x dt lC' ' The graphical interpretation of this elementary formula of calculus is illus trated in Figure [C.2. Substituting Equation (Cl) into (C.5) gives dt " {"u + dt} 37 (c'6) and along the curve defined by dx dt =u (C'7) the time rate of change of concentration is given by §£• = 0 (C.S) Equations (C.7) and (C.8) are the well-known Langrangian form of the advective transport equation and curve (C.7) obviously defines the three curves AB of Figure Cl. In deriving Equation (Cl) it is assumed that the variation of c and u with x and t is continuous. For the cases shown in Figure Cl, the vari ation of c is discontinuous and the partial derivatives 3c/3t and 3c/3x are not defined everywhere in the (x,t) plane. Thus, Equation (Cl) does not correctly represent the advection of a slug load (this is the fundamental reason for numerical dispersion). It is only when the solution to Equation (Cl) is ex amined along the curve (C.7) that the problem of the non-definition of the par tial derivatives is avoided. [See Equation (C.6)]. 142 if 3t = |f 8t Sx dt at ax dc _ ac + ac ^x dt at ax dt Figure C.2 Time Rate of Change of Concentration Along A Curve in the (x,t) Plane 143 Consider disturbing the system represented by Equation (Cl) by introducing an elemental slug load of Ac. This disturbance, or "information," will propagate through the system according to Equation (C.7). Thus, Equation (C.7) is of direct physical significance, and along this curve the solution to the advective transport equation reduces to the simple form of Equation (C.8). The fixed (Eulerian) space grid used to derive Equation (Cl) "masks" the under lying physics of the transport process and the simplicity of the solution. Now consider the equation of dispersive transport (C.2). If both E and A are constant in x and t, Equation (C.2) reduces to || - E 2i IC.9> ax The surface defined by Equation (C.9) for the dispersion of a slug load is shown in Figure C.3. The equation of this surface is [Fischer, 1966a] S(xft) ='-^~ exp{-'-—-} (CIO) 4IlEt where M is the mass per unit area;6f tracer released. At any timet1 the shape of the surface is a Gaussian distribution and c(x,t.') is uniquely defined by the location of the point of standard deviation. This is illustrated in Figure C.3. Thus, the locus of the standard deviation in the (x,t) plane directly des cribes the underlying mass transport process. From Equation (CIO) the locus of the standard deviation is the parabola x2 = 2E (Cll) or in the form of Equation (C.4) 144 Figure C.3 Dispersion of A Slug Load 145 If this system is disturbed by introducing an elemental slug load of Ac, the effect of this slug load, or its "information", will propagate through the system according to Equation (C.ll) or (C.12). Once again, the curve (C.12) is of fundamental significance and along this curve the concentration is given by _ M rl c(x,t) = -7=—-' (C.13) ^SHe X Finally, consider the hydrodynamic Equations (C.3). The well-known characteristics of this system are given by [Courant and Hilbert, 1962] §T = u ± (C.14) dt If this system is disturbed by introducing a small change in the water sur face elevation Ay, this "information" will propagate through the system accord ing to Equation (C.14) and along this curve Equations (C.3) reduce to [Henderson, 1966] d; dh . |u|u —- iu ± 2/gy} dt ---'w> = -9 37+J-1r <c-15> c y In summary, there are certain curves in the (x,t) plane that directly describe the propagation of information through the systems that Equations (C.l), (C.2) and (C.3) represent. These curves can be called charaeteristies of infor mation propagation, and along these curves the solutions to the respective equa tions are greatly simplified. This underlying simplicity is masked by the Eulerian nature of the derivation of the equations. The problems of numerical 146 dispersion and stability arise directly from viewing the numerical solu tions to Equations (Cl) to (C.3) from a fixed Eulerian space grid, rather than from a space grid along the characteristics. C.2 NUMERICAL DISPERSION Before discussing numericalddispersion in detail, it is necessary to briefly consider a few fundamentals of finite difference solutions to partial differential equations. In using finite difference methods to solve differential equations, derivatives are approximated by finite differ ence expressions and the differential equation is reduced to a difference equation. The justification of replacing the derivatives by difference ex^r-pressions can be illustrated by the Taylor's Series expansion of the variable U(x,t) with respect to t in the neighborhood of the point (x = jAx, t = nAt). This can be written as where the partial derivatives are assumed continuous and 0 < 0 < 1 ;[Richtmyer and Morton, 1967, p. 19]. By making At small enough, the difference expression of Equation (C.16) can be made to approximate the derivative to any desired degree of accuracy. The term on the right-hand side of Equation (C.16) is the truncation error that results when the derivative 8U/3t is replaced by the difference expression of Equation (C.16). When a partial differential equation in x and t is approximated by a partial difference equation, the truncation error consists of terms containing .n+1 (C.16) 147 Ax and At raised to various powers. The conditions under which the truncation error tends to zero is the problem of convergence. That is, as the mesh is refined (Ax, At -»• 0) , does the numerical solution converge to the true solu tion of the partial differential equation? The behaviour of thes truncation error as the solution progresses through time is the problem of stability. If the difference scheme is unstable, errors generated during the solution (for example, round-off errors) may become so magnified during the calcu lations as to make the final results meaningless. There are many fixed mesh difference schemes for approximating the derivatives of partial differential equations, but essentially they can be classified into explicit and implicit schemes. An explicit scheme uses forward time differences (see Figure C.4), and each difference equation con tains only one unknown variable which can be solved for explicitly. Implicit schemes use backward time differences (see Figure C.4) and each difference equation contains several unknown variables. To obtain a solution at any time step, the resulting system of simultaneous-iequations must be solved over the entire space grid. Generally, implicit schemes are unconditionally stable, whereas explicit schemes are at most conditionally stable. This is discussed further in Section C.3. Consider now the numerical solution of the one-dimensional mass transport equation (Cl). This equation is a component of the one-dimensional mass transport equation. Fixed mesh finite difference schemes generally do not simulate the advective transport process correctly, and result in an additional dispersive process being superimposed bn the actual advective and dispersive processes occurring in the river or estuary. This so-called numerical disper sion can be illustrated by considering the truncation error associated with 148 Forward Time Differences s Uin+I = function { u"-,, u" U" } -o-time = ( n + I )At time = n^t (j-l)Ax jAx (j+ I )AX Backward Time Differences : u n + I _ _ r n n*l n+l - function { Uj, Uj+|, U|_, } (J-I)AX —o— j AX (J+QI)AX time = (n + l)At time = nAt Figure C.4 Forward and Backward Time Differences 149 replacing the partial differential equation of advective transport (C.l) with its corresponding explicit difference equation ~n+l „n ,, »-,_n „n . n „n. C. - C. . (1-a (C . - C ) q (C - C .) where the flow is assumed to be uniform and steady, x = jAx, t = nAt and a is a weighting factor such that a = 0 for upstream differences, a = 1/2 for cen tral differences and a = 1 for downstream differences. These various differ ences are illustrated in Figure C.5. Eguation (C.17) is an example of an ex plicit difference scheme that uses forward time differences (the concentrations at time step n are used to determine the concentrations at time step n+1). A Taylor's Series expansion of the derivatives of Equation (C.l) shows that the difference equation (C.17) approximates the differential equation ar n ar n a2f n a2r n 2 2 O . + U(TT-) • = %UAx(l-2a) (±4) • " + O(Ax) + O(At)^ the terms on the right-hand side being the truncation error. As C is a solu tion of Equation (C.l) it also satisfies =0 (C.19) 9t 9x 2 Substituting Equation (C.19) into (C.18) and neglecting terms of order (Ax) . 2 (At) and higher, it is found that the difference equation approximates the differential equation 90 9C „ 9 C 9t * + Ep TI (C'20) 9x 150 Figure C.5 Upstream, Central and Downstream Space Differences 151 where E can be considered to be a coefficient of pseudo or numerical disper-P sion, and is given by E = %UAx{ (l-2a) -7^} (C.21) p Ax Thus, the solution to the difference equation (C.17) approximates the solution to an advective-dispersive equation rather than the correct ad vective equation. As Ax ->• 0, E^ •> 0 and the mass transport associated with the numerical dispersion becomes increasingly smaller.- However, in difference schemes Ax is always finite, and while the solution to difference equation converges to the solution of the correct advective equation in the limit (Ax = 0), it converges to the solution of an advective-dispersive equation at the level of application (Ax ^ 0). It is apparent from Equations (C.20) and (C.21) that numerical dis persion will always occur unless either U • 7^ = 1 - 2a (C.22) Ax or dX Ignoring the trivial latter case (the mass transport associated with any real dispersion is also zero), it is seen that numerical dispersion does not occur when upstream differences are used (a = 0) and Ax and At are defined by Ax TT = U (C.23) At 152 as was recognized by Bella and Dobbins [1968]. Under these conditions Equa tion (C.17) reduces to Cj+1 = Cj i = constant (C.24) Equations (C.23) and (C.24) are simply the finite difference versions of Equa tions (C.7) and (C.8). Thus, the use of upstream differences with the grid spacing defined by Equation (C.23) is equivalent to solving Equation (C.17) along the characteristic curve in the (x,t) plane. Under these conditions the underlying advective mass transport process is correctly simulated., and there is no numerical dispersion. Consider now using central differences (a = 1/2) to solve Equation (C.17). From-Equation (C.22) it is seen that numerical dispersion will always occur except in the trivial case of u = 0. According to Leendertse [1971b] there is no numerical dispersion when central differences are used. He apparently assumes that central differences correctly describe the advection process, and then uses this incorrect assumption to demonstrate that numerical dispersion will occur when upstream or downstream differences are used. The use of downstream differences (d = 1) will always result in numeri cal dispersion except when ff = -U (C.25) The use of downstream differences with Equation (C.25) defining the grid spacing can be interpreted as an attempt to determine upstream conditions from those downstream by working backwards through time (hence the negative sign) along the characteristic curve in the (x,t) plane. Although this procedure is im practical, it would correctly determine the preceeding concentration distributions 153 with no numerical dispersion if the final concentration distribution- was known. ': Bella and Grenney [1970] investigated the numerical dispersion result ing from the advection of a slug load by various fixed mesh finite difference;: schemes. They found that the coefficient of numerical dispersion was as given by Equation (C.21) and that no numerical dispersion occurred when up stream differences were used with Equation (C.23) defining the grid spacing. Numerical dispersion always occurred with central or downstream differences. Fox [1971] made a Fourier series analysis of the difference schemes of Bella and Grenney and showed that the effect of discretization was to introduce ampli tude and phase errors into the solution of the difference equations. His results demonstrate that neither amplitude nor phase errors occur when upstream differ ences are used with a grid spacing according to Equation (C.23). In conclusion, numerical dispersion occurs because fixed mesh finite difference schemes do not correctly simulate the advective transport process. (Although numerical dispersion has conly been demonstrated for explicit difference schemes, it also occurs in implicit difference schemes). The truncation error of the simple difference schemes discussed here corresponds to an actual physical mode of mass transport and will modify the advective transport and any dispersive transport that is occurring in a river or estuary. To correctly simulate the advective transport process, it is necessary to solve the advective mass transport equation along the characteristic curve in the (x,t) plane. The magnitude of the numerical dispersion depends directly 2 2 on the magnitude of 3 c/3x , and so is greatest for a slug load and less for a continuous release. A continuous release can be treated as a succession of 154 slug loads, and apparently the numerical dispersion due to any particular slug is compensated by the numerical dispersion of neighbouring slugs, as has been noted by Bella and Grenney [1970]. In deriving Equation (C.l) and writing the difference equation (C.17), it is implicitly assumed that c(x,t) varies in a continuous manner with x and t. This is certainly not the case for a slug load, as is seen in Figure C.l. It is only along the character istic curve that c(x,t) varies continuously, as is also apparent from Figure C.l. Thus, in using Equation (C.17) to advect a slug load, there is the added problem that the partial derivatives 8c/8t and 9c/9x are not defined, and the numerical dispersion can be interpreted as an attempt to define them. Numerical dispersion can be eliminated by solving the equation. advective-dispersive mass transport along the characteristics of advective information..propagatation. .Such a.^numerical solution is referred to as a characteristic method in Section 2.2, where the problem of the bookkeeping of solution results is discussed. Numerical dispersion can be controlled by the use of more sophisticated differencing schemes [see Fox, 1970], or by reducing the actual dispersion coefficients by E according to Equation (C.21). Numerical dispersion is apparently reduced in finite element solutions [Price et al.3 1968; Fox, 1970], but because of the fixed grid nature of finite element solutions, the numerical dispersion is probably not totally eliminated. C.3 STABILITY If a fixed grid difference equation is used to approximate a linear partial differential equation with constant coefficients, the von Neumann stability condition requires that the eigenvalues of the amplification matrix 155 should not exceed unity in absolute value for physically stable systems [Richtmyer and Morton, 1967]. To illustrate these ideas, consider the ex plicit finite difference scheme of Equation (C.17) that was used to approxi mate the advective transport equation. This difference scheme can be written cn+1 . c„ . _ «t,(1.al(c» . ^ + „(cn+i . cn)} (c26) where a = 0 for upstream differences, a = 1/2 for central differences and a = 1 for downstream differences. The von Neumann stability condition is based on a Fourier Series analysis of the difference equation. At any time, the concentration profile along the river or estuary can be considered to be composed of a number of Fourier components. The contribution of the k'th component to the concentration at point jAx at time nAt is given by „n+l „n r., .. i C . = C . exp i lki Ax i where i = /-l. Substituting into Equation (C.26) for the k'th Fourier component and simplifying gives n+1 „ n C = G C where G = l-0(l-2a)(1-cosO) + iBsinG (C.27) Ax and 0 = ikAx G is the amplification matrix of the difference equation (C.26). As C only depends on Cn, the eigenvalue of the amplification matrix is given by the value of G from Equation (C.27). Thus, the von Neumann stability condition 156 requires {l-8(l-2a)(1-cosO)}2 +{gsin 02} < 1 Evaluating this for the three cases of upstream, central and downstream differences, the stability criterion in each case is given by < Ax — or At > U ' (C If the stability condition is violated, certain Fourier components will be unacceptably amplified, and the value of the concentration, as deter mined by Equation (C.26) will oscillate with ever-increasing amplitude. Richtmyer and Morton [1967] give a very clear example of this oscillation for the simple dispersion equation (C.9). In the hydrodynamic equations, the instability manifests itself as an oscillation of the water depths and velocities (the water depths eventually becoming negative). Note that the von Neumann stability condition has been derived for linear equations with constant coefficients. This implies that the U of Equation (C.26) is con stant, or that the flow is uniform and steady. In a situation where U var ies with x and t, the stability condition is assumed to be determined by the worst case of Equation (C.28) along the estuary. Consider the simple dispersion equation (C.9). This can be approx mated by the explicit finite difference scheme cn+l _ cn = ^{c" - 2Cn + Cn > 3 3 (Ax)2 3+1 3 3-1 157 and if E is constant, the stability condition is given by [Richtmyer and Morton, 1967] 2EAt . < (Ax)2 T or -ll^— >_ 2E (C.29) The hydrodynamic equations (C.3) are non-linear and have non-constant coefficients. Richtmyer and Morton [1967] linearize the equations and obtain the following stability condition ^ > U± /qT At _ - <c-3°) It is apparent that Equations (C.28), (C.29) and (C.30) are the finite difference equivalents of Equations (C.7), (tc.ll) and (C.14), which describe the characteristic curve of information propagation in the (x,t) plane (see Section C.l). The physical reason for this can be seen from a consideration of Equation (C.26), the ^explicit finite difference approxima tion to the advective transport equation. This equation states that the value of C at the point jAx at time (n+l)At depends on the value of C at points (j-l)Ax, jAx and (j+l)Ax at time nAt. If this advective system is disturbed at point jAx at time nAt, as shown in Figure C.6, the disturbance propagates through the system according to Equation (C.7). If the stability condition of Equation (C.28) is violated, the disturbance will reach the point (j+l)Ax, 158 Stable Explicit Scheme uAt < Ax • Ac -o-Initial Position of Disturbance Position after At time = (n + I )At uAt U time = nAt Unstable Explicit Scheme uAt > Ax : Ac Position after At time = (n+ I )At Initial Position of Disturbance j-l uAt U Ax Ac time - nAt Figure C.6 Stable and Unstable Explicit Advective Schemes 159 and alter the concentration there, before the concentration at point jAx has been updated.; according to Equation (C.26) . In other words, the concen tration at point jAx at time (n+l)At now depends on what is occurring at a distance greater than Ax from jAx (see Figure C.6). Obviously, under these conditions the difference equation (C.26) is no longer physically meaningful, and it is only to be expected that it behaves in some strange manner. Simi lar reasoning holds for the relationship of Equation (C.29) to (C.ll) for the dispersive equation, and Equations (C.30) to (C.14) for the hydrodynamic equation. The problem of stability only arises with explicit differencing schemes, where the time differencing is forward. In implicit differencing schemes, where the time differencing is backward, the values are simultaneous ly determined at each grid point along the estuary. This technique automa tically adjusts for the effect of the information propagating distances greater than Ax in the time increment. A difference scheme is simply a?-means of determining the value CN+^~ from the value Cn. In effect, the difference equation maps C into C , and the eigenvalues of this mapping (in other words, the eigenvalues of the amplification matrix) define the characterise tic directions along which cn+^ is determined from Cn [see Sawyer, 1966]. It is only to be expected that these eigenvalues should reflect the underlying physics of the problem. C. 4 SUMMARY For any one-dimensional physical system, there are characteristic-" curves in the (x,t).plane that directly reflects the physical behaviour of the 160 system. These curves are often hidden or masked when a differential equa tion is derived relative to a fixed grid in the (x,'t) plane. The problem of numerical dispersion occurs because the equation of advective transport is solved relative to a fixed grid, rather than along its appropriate charac teristic. In an explicit difference scheme, the problem of stability arises for exactly the same reasons. If the equations are solved along their charac teristic curves (or by an implicit method) there are no stability require ments, and the relative size of Ax and At are determined by the variation in the dependent variable. APPENDIX D DETAILS OF THE SOLUTION SCHEMES OF THE HYDRODYNAMIC AND MASS TRANSPORT EQUATIONS D.l NUMERICAL SOLUTION OF THE HYDRODYNAMIC EQUATIONS The fixed mesh, explicit finite difference method of Dronkers [1969] was used to obtain a numerical solution to the hydrodynamic eguations. The finite difference forms of Equations (3.1) and (3.2) are respectively 2n+l 2n-l .At , . 2n+l,..2n-l 2n-l, At,,2n .2n -» U2m = U2m " (2Ax"} U2m (U2m+2 " U2m-2) ~ gAx-(h2m+l ~ h2m-i) cV" (D.l) Y2m and 2n+2 _ 2n At 2n+l 2n 2n+l 2n n2m+l " h2m+l " . _2n (U2m+2 A2m+2 " U2m' A2m) (D,2) ^xb^ 2m+l where and h^ are respectively the mean velocity and surface elevation at the grid point given by x = mAx at time t = nAt. Note that the difference scheme employs central differences and that the "bar" sign has been dropped from the cross-sectionally averaged variables and parameters to avoid confusion with the subscripts and superscripts'. (This convention will be followed when presenting finite difference quantities). Because an explicit finite difference scheme was used to solve the hydrodynamic equations, the relative size of Ax and At is governed by the 161 162 stability criterion Ax • -• • r-=y At-'U ± /gY This stability criterion is discussed in detail in Appendix C. Because of the friction term, the hydrodynamic equations are naturally dissipative, and this helps maintain stability. In fact, stability was found to be de pendent on a minimum value of friction. For the cross-sectional gemoetries and flow and tidal conditions of Section 4.1.4, the hydrodynamic equations became unstable when Manning's "n" was less than 0.012. (Manning's "n" was constant throughout the estuary for this investigation). The stability cri terion of Appendix C is independent of the effects of friction. However, it was derived from a linear stability analysis on the non-linear hydrodyna mic equations. The true non-linear stability criterion may be more strict than the derived linear criterion, and this is reflected by the necessity of a minimum level of friction to preserve stability. The solution points for the finite difference scheme are shown in Figure D.l. It is seen that the velocity and elevation points are stag gered in both time and space. In effect, the values of velocity and water surface elevations are marched forward through time in a "leap-frog" manner. From Equation D.l, the velocities at any time step 2n+l are determined by the surface elevations at time step 2n and the velocities at time step 2n-l. The surface elevations at time step 2n+2 are then determined from the velo cities at time step 2n+l by'using Equation•(D.2). The fixed mesh of space points or ";stations" used in solving the hydrodynamic equations is shown in Figure B.4. Velocities are evaluated at odd -numbered stations and water surface elevations are evaluated at even-2n+2 2n+l 2n 2n-l u 1 u s A u 1 % t ' > 1 h X ,h ? u A V A / U \ i u % f % h ? h f \ ? U (> < > u c u > 2m-3 2m-2 2m-| 2m 2m+| 2m + 2 x longitudinal distance Figure D.l Explicit Finite Difference Grid of the Hydrodynamic Equations [After Dronkers, 1969] 164 numbered stations. At channel junctions, it is necessary that the junction station be a water surface elevation station, and this accounts for the non-coincidence of the junctions of the model estuary and real estuary in several cases (see Figure B.4). In applying the hydrodynamic model, a design tidal cycle of selected flow and tidal conditions is chosen. The freshwater inflow forms a boundary condition at Chilliwack and the tidal rise and fall of the water surface forms boundary conditions at the seaward ends of the four channels that emerge from the delta. The velocity and water surface elevations through out the estuary are.-rset to initial values, and Equations >(D.l) and (D.2) are then used to march the velocities and water surface elevations through time. The initial values of velocity and water surface elevations do not have to be exact as the model will converge to the true initial values after several tidal cycles.; The hydrodynamic model was programmed in high speed FORTRAN (FORTRANH) for solution by digital computer. The program consisted of 670 active state ments and required approximately 40 seconds to analyse two complete tidal cycles. (The first tidal cycle was required for the estuary to converge to the true initial values of velocity and water surface elevation). The values of velocity and cross-sectional area were recorded at half-hourly intervals and later used in solving the tidally varying mass transport equation. D.2 NUMERICAL SOLUTION OF THE TIDALLY VARYING MASS TRANSPORT EQUATION From Section 3.2 it is seen that the solution of the tidally vary ing mass transport equation over any time increment involves an advection step, a dispersion step and a source-sink step. The finite difference form 165 of Equation (3.4) used to advect the moving points along the characteristics during the advection step was n+1 n n.. x_. = x.. + u At (D.3) where n+1 . x. is the position of moving point j at the end of time -1 increment n; x" is the position of moving point j at the start of 2 time increment n; and un ' is the average velocity between positions x1? and x1?"*"^" during time increment n. ^ ^ The hydrodynamic model was used to obtain the longitudinal velocities and cross-sectional areas throughout the estuary at half-hourly intervals during the tidal cycle. The value of un can be obtained, from these velocities. In the dispersion step, the concentration of the moving points is adjusted for the effects of dispersion during the time increment. Both an explicit and an implicit finite difference scheme were investigated for the dispersive step. The implicit scheme was the Crank-Nicholson scheme described in Richtmyer and Morton [1967]. Irrespective of whether an impli cit or explicit scheme is used, the "information" propagates through a dis persive system according to x2 = 2Et (D.4) as is discussed in Appendix C. The finite difference form of Equation (D.4), namely *< <§>2 defines the stability criterion of an explicit scheme/ as is also discussed in Appendix C. Equation (D.5) also determines the response of an implicit system, and although the implicit system is unconditionally stable, the con vergence criterion is related to Equation (D.5). In other words, if At is significantly larger than the At of Equation (D.5) an implicit scheme may converge to the wrong solution. Of the two schemes, the implicit scheme was slightly faster, but because of the large somewhat ill-conditioned matrices involved (of order 150) there were uncertainties in the signifi cance of round-off errors. Consequently, the simpler explicit scheme was used. (The ill-conditioned nature of these matrices was due to the variable spacing of the moving points. This variable spacing is discussed in Chapter The following explicit central difference equation was used to Equation (3.6) 3). n+1 n c. + 2At 3 A.(Ax) . 3 3 n+1 j+l,j-l Ax'j,j-l^j (D where c. is the concentration of moving point j at the start of ^ time increment n; j,j-l n = x. - x 3 n j-l x. is the position of moving point j at the start of time 3 increment n; (EA)n j,j-l is the average value of EA between moving points j and j-l during time increment n; xj,'[ and A. is the average cross-sectional area at moving point j during time increment n. 167 Equation (D.6) is the usual explicit central difference approximation except that it is applied over a grid where Ax is not constant. The stability re quirement of Equation (D.5) generally results in a At smaller than the basic time increment of one hour. When this occurred, Equation (D.6) was solved repeatedly within the hour for as many iterations required by Equa tion (D.5). (This assumes that the relative spacing of the particles and the values of E and A remain constant during the hour — which is a reason able assumption). Finally, in the source-sink step, the concentrations of the moving points are adjusted for the effects of any source-sink processes occurring. The following finite difference equation can be used to approximate the source-sink equation (3.7) +1 n c. = cP + At .Z.S. (D.7) 3 3 1=1 i although it is noted that Equation (3.7) can be solved analytically during the time increment. The tidally varying mass transport model was programmed in high speed FORTRAN (FORTRANH) for solution by digital computer. The program con sisted of some 1,300 active statements and required approximately 100 seconds to analyse six tidal cycles. D.3 NUMERICAL SOLUTION OF THE TIDALLY AVERAGED MASS TRANSPORT EQUATION The one-dimensional, tidally averaged mass transport was solved by the method of Thomann [1963]. In this solution, the estuary is divided into a number of segments or boxes, as shown in Figure D.2, and each segment is assumed to be completely mixed. If the tidally averaged transport processes 169 and waste discharges are steady in time, a inass balance over segment i for a substance undergoing firsfcrorder decay gives (see Thomann [1971] for de tails) Q{a. , .c, + (1-a. . .)c.} - Q{a. _.c. + (1-a. . Jc .} * i-l,i i-1 i-l,i i i,i+l i i,i+l i+l + E. .(c. -c.) + E. . ,(c. -c.) - K.c.V. + W. = 0 (D.8) i-l,i i-1 l i,i+l i+l l ill l where c is the concentration in segment i; V^ is the volume of segment i; W. is the mass of waste substance discharged into segment i per tidal cycle; K^ is the decay coefficient for segment i; Q is the tidally averaged discharge through the estuary (the freshwater discharge); a. is the tidal exchange coefficient between segments i 1,1+1 and (i+l); and • E^ is the "effective dispersive" transport between segments lf i and (i+l). The subscript notation of the various terms is illustrated in Figure D.2. The first two terms on the right-hand side of Equation (D.8) are the tidally averaged advective transport into and out of segment i. The factor a is a weighting factor used to determine the concentration at the interface of two segments from the concentration at within each segment. In purely tidal flows a is set equal to 0.5 to allow for the effects of flow reversal, whereas in a river flow situation a is set equal to 1.0 as the flow is always downstream. The next two terms on the right-hand side of the equation 170 represent the net dispersive transport of mass into segment i from the neighbouring segments. E' is given by • "i,i+l E- _ Ej,i+1 Aj,i+1 where and Li,i+1 E^ is the effective coefficient of dispersion over a ' tidal period at the interface of segments i and (i+1); A. is the cross-sectional area (tidally averaged) of the ' interface between segments i and (i+1); is the average of the lengths of segments i and Xi+D• The final two terms on the right-handclside of Equation (D.8) represent the effects of decay and waste discharge. An equation similar to (D.8) can be written for each of the n segments of the estuary to give a system of n simultaneous, linear, difference equations. Equation (D.8) can be written where a. . .c. , + a..c. + a. . .c.^. = W. (D.9) i,i-l i-l ii I i,i+l i+I l a. . .. = -a. . . Q - E. , . ; i,i-l i-l,i i-l,i a-.... = Q{a. . - (1- c. , .} + E. + E. + V.K.; ii i,i+l i-l,i i-l,i i,i+l i l a. . . = (1-a. )Q - E. i,i+l i,i+l i,i+l. In matrix notation, the system of Equations (D.9) can be written k £ =  (D- 10) where A is a (nxn) tri-diagonal matrix and C and W are (nxl) column matrices. ri, • J O, r\j ' 171 Thomann [1971] gives details of the complete BOD-DO system of equations. Thomann's model is essentially a finite fixed grid finite dif ference model that uses central differences for the dispersive transport, and central differences for the advective transport when a = 0.5 and up stream differences for the advective transport when a = 1.0. (Upstream and central differences are described in Appendix C). Thomann's model is simi~-; lar to an implicit finite difference scheme. In both schemes the concentra tion at a gridpoint or in a segment depends on the concentrations at neigh bouring gridpoints or segments, and consequently the response of the estuary is determined by a square tri-diagonal matrix in both equations. Finite element solutions are also governed by a square tri-diagonal matrix, as is discussed in Section 2.2.3. Thomann's difference equations are unconditionally stable and do not suffer from numerical dispersion, as is discussed in Chapter 3. How ever, there is a non-negativity requirement for each segment given by E. a.,i+l > 1 - X' (D.ll) i Q If this criterion is violated the discharge of a waste substance into seg ment i results in a negative concentration in the segment. The physical reason for this is that more substance is being transported out of the segment per tide cycle than is being added. The substance is transported upstream by dispersion and downstream by advection and dispersion. Rearrang ing Equation (D.ll) gives E. . n > (1-a. . ,,)Q i,i+l i,i+l x which imposes limits on the relative size of the dispersive and advective transport processes. 172 Finally, the ease with which Thomann's approach handles the separate channels of the estuary should be mentioned. Figure D.3 shows the matrix A of Equation (D.10). Note that the three channels of the estuary are contained in the one matrix. In effect, the matrix is partitioned into three separate blocks, each of which represents a single channel. Note that each block or channel is uncoupled from the others except at the junction stations, where an additional term, which is not tri-diagonal, appears in the rows and columns of the matrix. These additional terms reflect the extra boundary through which mass trans port occurs at the junctions. The tidally averaged mass transport model was programmed in FORTRAN for solution by digital computer. The program consisted of some 250 active statements and required approximately 10 seconds to determine the steady state response of the estuary. 173 Figure D.3 The Matrix A of Eguation (D.8) For The Fraser River Estuary APPENDIX E ESTIMATION OF LATERAL DISPERSION Existing theories of lateral dispersion were used to estimate the time of cross-sectional mixing in the Main Arm - Main Stem of the Fraser River Estuary. The predicted time of cross-sectional mixing appears high for the conditions of this study. More recent work on secondary currents in rivers was used to develop revised estimates of the cross-sectional mixing time. These revised estimates indicate much faster mixing over the cross-section. E.l EXISTING ESTIMATES OF LATERAL MIXING For steady uniform flow in straight channels, the coefficient of lateral dispersion is given by [Fischer, 1969a]. ez = 0.23yU^ (E.l) where y is the mean depth of cross-section; and is the shear velocity. The presence of bends in the channel of a river or estuary induces spiral secondary currents that increase the rate of lateral mixing. For steady flow around a bend, the increase in the coefficient of lateral dispersion due to these secondary currents has been estimated as [Fischer, 1969a] -2-3 Se = - U/ • I (E.2) z '5 2 k R U* 174 I 175 where u is the mean longitudinal velocity; k is Von Kantian's constant; R is the radius of curvature of the bend; and I is a factor (negative) that depends on channel friction and k, and is evaluated by Fischer. (A typical value of I is -0.3). Ward [1972] measured the coefficient of lateral dispersion from laboratory experiments of oscillatory flow (purely tidal) in a channel con sisting of a series-of bends. His results are of the form e' = ayU* (E.3) z * where and e is the tidally averaged coefficient of lateral dispersion due to oscillatory flow; is the tidally averaged shear velocity due to oscillatory flow; a is a factor depending on the ratios of the depth and width of the channels to the radius of curvature of the bends. (For straight channels a = 0.5, and for wide shallow channels with a small radius of curvature a^2.0). In a tidally varying situation where there is both a steady velocity compo nent and a sinusoidally varying component, the tidally averaged shear velocity can be estimated from where h is Manning's "n"f is the steady velocity component; 176 and U is the amplitude of the oscillatory velocity component. Equation (E.4) is obtained from the Manning formula relating steady velocity to frictional effects. According to Odd [1971], the Manning or Chezy formu lation should be an adequate description of frictional effects in fast flow ing, well-mixed estuaries. The following estimation of dispersion coefficients is limited to the Main Arm - Main Stem of the Fraser River Estuary. This is the widest channel and carries the bulk of the flow through the delta. Along this chan nel, there are 12 significant bends whose radii of curvature range from 7,000 feet to 35,000 feet, as illustrated in Figure E.l. The variation of longi tudinal velocity throughout the tidal cycle at three stations along the Main Arm - Main Stem is shown in Figure E.2. The predominance of the tidal component in the lower reaches and the steady component in the upper reaches is apparent. The following procedure was used to estimate the coefficients of lateral dispersion. The velocity at each bend was divided into a steady com ponent Uf, which is given by the freshwater discharge through the tidally averaged area, and an approximating sinusoidal component of amplitude Ut. Equations ;(E-;.i) and (E.2) were used to estimate the coefficient of lateral dispersion due to the steady component of velocity ('e ) , and Equation ;[E.3] was used to estimate the coefficient of lateral dispersion due to the sinu soidal or oscillatory component (e^). Equation (E.4) was used to estimate f t the shear velocities due to the steady (U^) and oscillatory (U^) velocity components. The total lateral dispersion was assumed to be the sum of these two separate effects. The results of the calculations are shown in Figure E.l Bends Along The Main Arm - Main Stem Q= 36,500cfs at Chilliwack Tidal Range at Steveston Jan.24,1952 178 24( Hours) 24( Hours) Stat.No.33 24(Hours) 24(Hours) -2l Figure E.2 Tidally Varying Velocities TABLE E.l ESTIMATION OF COEFFICIENTS OF LATERAL DISPERSION BEND NO. STAT. NUMBERS +J <D £ 0) fa -p fa o 0) w \ iw • ID -P fa o <D W \ D fa ~G • V <D Ul \ <4-l * • D -P fa « o <u w \ +1 * • D -P fa o <u w fa \ A \ o •P * D 1 >i \ -P (J M-l * P 1 >i \ u u CO W \ fa 1 57-60 17,000 24 1.2 0.1 0.032 0.088 0.005 0.088 0.12 0.14 1.0 0.24 0.6 2 55r57 7,000 24 1.2 0.8 0.032 0.088 0.037 0.097 0.29 0.34 1.9 0.29 2.3 3 52-55 7,000 35 0.9 0.7 0.032 0.062 0.031 0.070 0.29 0.40 1.8 0.37 2.8 4 46-52 17,000 31 0.9 1.2 0.032 0.063 0.054 0.087 0.13 0.18 1.0 0.24 2.2 5 41-46 16,000 37 0.8 1.3 0.032 0.055 0.057 0.083 0.13 0.23 1.0 0.26 2.7 6 36-41 35,000 32 0.7 1.3 0.032 0.049 0.058 0.081 0.06 0.09 0.6 0.23 1.5 7 32-35 7,000 43 0.7 1.5 0.025 0.036 0.050 0.066 0.29 0.62 1.8 0.62 4.8 8 29-31 7,000 38 0.7 2.2 0.025 0.037 0.074 0.091 0.29 0.54 1.8 0.52 5.8 9 19-23 16,000 38 0.6 2.5 0.025 0.033 0.088 0.10 0.13 0.24 1.0 0.28 4.7 10 14-17 8,000 38 0.7 3.3 0.025 0.037 0.11 0.12 0.25 0.48 1.6 0.45 7.3 11 8-14 33,000 34 0.6 3.0 0.025 0.032 0.10 0.11 0.06 0.10 0.6 0.24 2.3 12 2-8 27,000 34 0.6 3.5 0.025 0.032 0.12 0.13 0.07 0.13 0.6 0.24 2.7 180 Table E.l, the combined coefficient of lateral dispersion being designated e . c From the results of Table E.l, the average value of the coeffi cient of lateral dispersion along the Main 'Arm - Main Stem is approximately 3.3 square feet per second. Assuming an average width and depth of channel of 1800 feet and 30fcfeet respectively, and an effluent discharge at the channel edge, the time required for 80 per cent mixing in the lateral direc tion is 55 hours [Ward, 1973]. (The percentage lateral mixing is defined as the ratio of the lateral root-mean-square concentration deviation to the aver age lateral concentration). This estimate of 55 hours for the time of 80 per cent cross-sectional mixing seems very high. Strong secondary currents around bends 9, 10 and 11 are observable in the estuary during ebb flow conditions. In fact, if the fishermen lose a net in the Main Stem, it is often washed ashore on the north bank of bend 11 by these secondary currents. In a surface float study during ebb tide conditions, the floats were also washed ashore on the north bank of bend 11. It seems that the results of Table E.l underestimate the effects of the secondary currents on the lateral mixing process. In view of this, an attempt was made to estimate the velocities of the secondary currents and the influence of the secondary currents on lateral mixing. E.2 VORTICITY ESTIMATE OF SECONDARY CURRENTS First consider the lateral mixing due to secondary currents. For the sake of simplicity, the vertical distribution of the secondary velocities is assumed to be linear, as shown in Figure E.3. This velocity distribution should be a satisfactory approximation for wide, shallow cross-sections. 181 w(?7) = Ws( !- 2 77)1 V = y/y _ Ws  ^^^^^^^^^^^^ B Figure E.3 Assumed Linear Distribution of Secondary Velocities 182 The variation of velocity with depth is given by w(n) = Ws(l - 2n) (E.5) n = y/y where Wg is the surface velocity of the secondary current. In the manner of Elder [1959] and Fischer [1969a]the coefficient of lateral dispersion can be estimated from e = -yr2J1w(n){fT,w(Ti> [Jn — dn]dn>dn (E.6) z . ' e o o o y where e^ is the coefficient of vertical mixing. If we assume the flow to be steady and. uniform, the vertical profile of the longitudinal velocity u is logarithmic, and is given by [Einstein, 1972] -, , U* • r9.05yU*,, u(y) = r- ln{ £ }} (E.7) is. V where y is measured vertically upwards from the bend; and v is the kinematic viscosity. For this logarithmic profile, the value of is given by [Fischer, 1969a] ey =: k (l-n)ny<U* (E.8) where k is von Kantian's constant. Using the relations of Equations (E.5) and (E.8), the integration of Equation (E.6) gives e w 2 -2— = 0.42{-S-} (E.9) y U* 183 This result closely agrees with the value of e W -2— = 0.5{-^}2 y U* U* that Ward [1972] gives for the same distribution of secondary velocities. Thus, it only remains to estimate or measure the secondary velocity Wg to determine the coefficient of lateral mixing. Quick [1973] has developed a theory of river meanders that nicely predicts the types of meanders that are observed in nature. The existance of secondary flows when bends or meanders are present in a river is well known, and Quick's theory has been used to obtain an estimate of the secondary velocity Wg. The basis of this theory is the generation of vorticity by the gradient of the vertical distribution of longitudinal velocity u, and the sub sequent advection of this vorticity by u. This is illustrated in Figure E.4, the sense of the vorticity being given by the right-hand-screw rule. Because of the large velocity gradients of u close to the bed, the vorticity is gener ated in this region as intense, small-scale, highly anisotropic vortices. The vorticity then diffuses vertically into the flow and is advected along by u. During this diffusion and advection, the vorticity degenerates into increasing ly isotropic turbulence and is ultimately dissipated by frictional effects. However, in any steady flow situation there must be a balance between the generation and dissipation of vorticity. The vertical logarithmic profile of u reflects the energy balance of a steady flow situation, and presumably the vorticity balance also. Because of the velocity gradients of this pro file, the flow at any section has a "bulk" vorticity given by P = 9u(y) 3y (E.10) average 184 Figure E.4 Interaction of u and £ In A Straight River 185 where the velocity gradient is averaged over the depth of flow. If this argument is correct, there is an intermediate stage between the generation and dissipation of the vorticity, in which the vorticity retains an organized, anisotropic state. This state is reflected by the velocity gradients of the logarithmic velocity profile. In a straight river, the velocity vector U and the vorticity vec tor are perpendicular to each other and do not interact. This is illus trated in Figure E.4. However, the situation is different when the flow moves around a bend. Initially, the water on the inside of the bend flows faster than the water on the outside of the bend, but in the bend proper, the faster filaments of flow move to the outside of the bend, and the outside water moves faster than- the water on the inside of the bend. This is illus trated in Figure E.5. At A the faster filament is on the inside of the bend while at B it is on the outside of the bend. Because of the lateral velocity gradient of u as the flow moves around the bend, the vorticity is being gener ated and advected at an angle 0 to uy as illustrated in Figure E.5. Note that 0 varies as the flow moves around the bend and this variation of 0 will depend on the radius of the curvature of the bend. Quick's theory predicts this angle?© between u and £ . Thus, there is a streamwise component of vorticity given by E, = k |cos0 (E.U) X z and it is this streamwise component that gives rise to the secondary current as the flow moves around the':>berid. Substituting the logarithmic profile (E.7) into (E.10) and aver aging over,.the.;depth.between the limits y = 11.0v/U* and y =,y, the lower limit being an estimate of the depth of the laminar sublayer, gives 186 Figure E.5 Interaction of u and 5 Around A Bend 187 aU* K = (E.12) where a = ln{r-T-—r} (E.1311. Ov and thus from Equation (E.ll) = aUJ*cos0 ;(E;14) k y Now, consider the circulation around the section ABCD of Figure E.3 due to this streamwise vorticity. This circulation is given by «r = <j> ? • dA X ABCD ^X = i- d U*b cosG (E.15) where b is the width of the river. Now, the circulation around the cross-section ABCD is also given by r = <j> U • dl (E.16) X ABCD ^ ^ Assuming there is sufficient time for the streamwise vorticity to "arrange" itself into the linear distribution of secondary velocities of Figure E.3, Equation (E.16) can be evaluated as T ~ 2W b (E.17) x s where the contributions over BC and DA have been ignored (the sectidn is assumed to be wide and shallow). Finally, putting k = 0.4 and substituting Equation (E.15) into (E.17), the following expression for Wg is obtained W = 1.25aU. cosG (E.18) s * t 188 The results of a surface float study indicate a cross-channel e surface velocity of approximately 0.6 feet per second around bend No. 11 of Figure E.l. The study was made on March 20, 1973 when the freshwater discharge at Hope was 31,000 cubic feet per second. The floats had drogues to minimize the effects of wind and were released during the initial period of the strong ebb phase of the tide. The following velocity values were ob tained : u =3.1 feet per second (from floats); and = 0.15 feet per second (Manning's equation) In his study of stable meanders, Quick found that for geometrically similar meanders to bend No. 11, the value of 9 varied from about 50 to 70 degrees. Taking as above and y = 34 feet (Table E.l) -5 w = 1.6 x 10 square feet per second; and 0 = 60 degrees we obtain from Equations (E.13) and (E.18) Wg ? 0.7 feet per second which agrees well with the value obtained from the float studies. Assuming a linear variation of secondary flow with depth (as in Figure E.3) and using a value of Wg of 0.6 feet per second, the coefficient lateral dispersion from Equation (E.9) is given by 189 e — = 6 (E.19) ? U* or ez = 30 square feet per second. The value of Equation (E.19) is considerably larger than other : values that have been reported (for example, the values in Table 1 of Ward [1972]), and it should be emphasized that the estimated value has not been confirmed experimentally. A possible reason for this high value is that the Fraser is essentially a "tidal river" with both High -freshwater and tidal flows. The Fraser is relatively wide when compared to other rivers where values of ez/y U^ have been measured, but relatively narrow when com pared to estuaries (see Table 1 of Ward). It is noted that the dispersion equation (E.9) for determining the lateral dispersion coefficient in second ary flows is very sensitive to the value of (Ws/U^). Using a value of = 30 square feet per second, the previous estimate of 55 hours for 80 per cent cross-sectional mixing is reduced to five hours. This will underestimate the time of cross-sectional mixing as the flow is not steady during the tidal cycle. Because of the highly assymetrical nature of the tides, there is only one strong ebb and flood tide in each double tidal cycle of 25 hours to generate secondary currents (see Appendix B for typical tides). Thus, within any tidal cycle, the later al mixing will vary from a maximum during the strong flood or ebb to a minimum at times of slackwater. Consequently, an estimate of the effect of cross-sectional mixing is probably one tidal cycle (12.5 hours) for the lower reaches of the estuary and Pitt River and 1-2 tidal cycles for the upper 190 reaches where the tidal effects are smaller. On this basis, the effective or tidally averaged coefficient of lateral dispersion would be about 15 square feet per second in the lower reaches and seven square feet per second in the upper reaches of the Main Arm - Main Stem. On the basis of the relative magnitudes of the tidally averaged values of y and U^, the lateral dispersion in the North Arm is estimated to be five square feet per second and 10 square feet per second in Pitt River. At this stage, the theory relating vorticity and secondary currents is neither fully developed theoretically nor confirmed experimen tally, but future work is planned in both directions. In concluding, it is noted that some values quoted for lateral dispersion coefficients are based onithe results of laboratory experiments. Cross-rsectional mixing in the presence of secondary currents is a complex three-dimensional pheno mena, and it may be that some components of this process are not being scaled properly in model experiments. APPENDIX P ESTIMATION OF LONGITUDINAL DISPERSION The coefficients of longitudinal dispersion due to the effects of vertical and lateral velocity gradients are estimated for the Fraser River Estuary. Simple approximations are given for the time-dependent behaviour of the longitudinal dispersion coefficient during the initial period before cross-sectional mixing is complete and during the tidal cycle. The predicted tidally varying concentrations during the first double tidal cycle were found to be very sensitive to assumptions about the time-dependent behaviour of the dispersion coefficient. Neither the magnitude nor time-dependent behaviour of the dispersion coefficients has been verified by field measurements, and it is recognized that they may be in error. F.1 GENERAL When effluent is discharged into a river or estuary it is dispersed in the longitudinal direction by the effects of both vertical and lateral velocity gradients, as described in Appendix A. The effluent must disperse over a reasonable depth and width of the estuary before these velocity grad ients can exert a significant effect on the longitudinal dispersion process. Most rivers and estuaries are much wider than they are deep and consequently mixing in the vertical direction is much more rapid than in mixing in the lateral direction. Thus, when a "parcel" of effluent is initially discharged, 191 192 the longitudinal dispersion is essentially due to the effects of vertical velocity gradients. However, as the effluent spreads across the cross-section, the lateral velocity gradients exert an increasing influence on the longitudinal dispersion process. In an estuary, the longitudinal dispersion is complicated by the effects of tidal flow reversal;,. If the cross-sectional mixing is not essentially complete within a tidal cycle, some of the longitudinal dispersion due to the oscillatory flow will be "undone" by the effects of flow reversal [Holley et al.3 1970]. The Fraser River Estuary falls into a class that Holley et al. [1970] describe as "sinuous, multi-channeled or island-studded estuaries." Because of the complicating effects of the bends, islands and junctions of the estuary, the coefficients of longitudinal dispersion can only be reliab ly determined by field dye studies. Time and expense precluded such studies, and in the absence of adequate field data, the work of Fischer [1966b, 1969b] and Holley et al. [1970] has been used to obtain preliminary estimates of the longitudinal dispersion. F.2 TIDALLY AVERAGED COEFFICIENTS OF LONGITUDINAL DISPERSION For steady flow the longitudinal dispersion due to the effects of vertical and lateral velocity gradients can be respectively estimated as E -2— = 6 (F.l) y u* 193 where u' is the square of the velocity deviation and is given by 2 - 2 u' = (u(y,z) - u) the over-bar in Equation (F.2) signifying the average cross-sectional value. Eguation (F.l) was obtained by Elder [1959] for the logarithmic velocity profile. In obtaining Equation (F.2), Fischer : [1966b] assumed the lateral mixing to be due to turbulent diffusion alone. According to the triple integral of Equation (E.6), the longitudinal dispersion due to lateral velocity gradients will vary inversely as the coefficient of later al dispersion. This is because mixing over the cross-section tends to "undo" the effects of dispersion in the longitudinal direction. The lateral dis persion in the Fraser River estuary is thought to be quite high because of the effects of secondary flows (see Appendix E), and Equation (F.2) is sub sequently adjusted for this effect. To apply Equation (F.2) it is necessary to have some estimate 2 of u' and . Two velocity profiles were made in the Main Arm at Stations Nos. 14 and 15 on April 4, 1973. The freshwater flow at Hope was 34,700 cubic feet per second, the tidal range at Steveston was 11 feet and the measurements were taken during the strong ebb phase of the tide. The depth-averaged velocity profiles across the sections are shown in Figure F.l, and for these conditions, the various velocities are estimated to be u = 3.6 feet per second (measured); = 0.18 feet per second (Manning's equation); = 0.9 feet per second; and U. = 3.3 feet per second Q = 31,000 cfs (Hope) April 4,1973 194 500 1000 Distance from North Bank ( feet ) 500 1000 1500 Distance from North Bank (feet) Figure F.l Lateral Velocity Profiles At Stations No. 14 and 15, Main Arm 195 where it is assumed that the total velocity u consists of a steady compo nent (freshwater) of magnitude Uf and a sinusoidally varying oscillatory component.(tidal) of amplitude Ufc. Stations Nos. 14 and 15 are located in bend No. 10 of Figure E.l, a region of strong secondary flows. Using the average of the measured -2 values of u' at both stations, Equation (F.2) can be evaluated as E — = 3600 (F.3) which is much higher than recorded values in other rivers [Fischer, 1966b, Table 1]. However, the value of Equation (F.3) has to be reduced to account for the greater cross-sectional mixing in the Fraser. The relative magnitude of lateral mixing due to turbulent diffusion is 0.23, as in Equation (E.l), whereas the relative magnitude for the effects of secondary currents in the region of the estuary has been estimated to be 6.0 (see Appendix E). Thus the adjusted value of Equation (F.3) is E —— = 140 (F.4) y which is much more reasonable when compared to the values that Fischer lists. Equation (F.4) is an estimate of the "instantaneous" longitudinal dispersion due to the lateral velocity gradients existing at one particular time during the tidal cycle. According to Fischer [1969a], the effects of the steady and oscillatory components on the longitudinal dispersion are 196 separate and additive, and thus Equation (F.4) can be separated into the following steady and oscillatory components 5"f = 140 (F.5) and = 140 (F.6) where For the sake of simplicity, the subscript z has been dropped. Note that Equations (F.5) and (F.6) sum to give the correct combined dispersion of Eguation (F.4). To apply the analysis of Holley et at. [1970], the oscilla tory component has to be corrected back to a tidally averaged value. In the manner of Equation (E.4), this is estimated to be Et ^ = 120 (F.7) y & * In the absence of other field data, Equations (F.5) and (F.7) have been used to estimate the longitudinal dispersion due to lateral velocity gradients throughout the estuary. The equations are applied as is for the Main Arm - Main Stem of the estuary, but have been corrected for the different depths and widths of the other channels to give E_ - f y u* = 75 E } North Arm (F.8) = ,65 y^ 197 = 120 Pitt River (F.9) Y< Note that there is no steady dispersion component in Pitt River, the flow is purely oscillatory. The estimates of the various velocity components and the steady and oscillatory dispersion components are given in Table F.l. If the cross-sectional mixing is not essentially complete within a tidal; cycle, some of the oscillatory dispersion is undone by the effects of tidal flow reversal. To account for this effect, the values of E^ have been reduced to their effective values E£ according to the procedure in Holley et at. These results are also shown in Table F.l, the effective per iod of flow oscillation is designated T and the ratio of this value to the time scale of lateral mixing is designated T'. Finally, the steady and effective oscillatory dispersion components are summed to give a combined tidally averaged coefficient of longitudinal dispersion due to the effects of lateral velocity gradients (Ec). It is recognized that the values of the dispersion coefficients given in Table F.l have been obtained from very limited field data and de pend heavily on the unverified estimate of lateral dispersion from Appendix E. Consequently, the estimates may be considerably in error. However, the ratio Ef/y and the absolute value Ef are reasonable when compared to values measured in other rivers [Fischer, 1966b, Table 1] and the maximum value of the effective oscillatory component Efc of 450 square feet per second agrees well with the upper limit of approximately 500 square feet per second that Fischer [1969b] suggests. The combined dispersion coefficient Ec shows an TABLE F.l ESTIMATED COEFFICIENTS OF LONGITUDINAL DISPERSION DUE TO LATERAL VELOCITY GRADIENTS (January 24, 1952: Freshwater Discharge at Chilliwack 36,500 cfs; Tidal Range at Steveston 11 feet) . . u u • . - . . o o CJ CU CU u * o o o o cj CU cu — CO co D •H cu cu 'CO cu CO co co CO \ \ w • co co 10 CO co \ \ \ U CM CN 1 >i z to \ \ \ \ \ CN CN CN - 4J • u • \ z • 4-> • M-l * • 4-> * • o * • IW • 4J • N • o - W IP W 4J O < D -P D 4J D -P , D 4-> P 4J W 4J fa V Ul, 4-> En a EH fa fa w u EH CO fa ' fa fai fa fa_ fa fa 1-10 0.6 3.5 0.025 0.11 0.12 120 450 15 25 1.7 450 570 140 11-18 0.6 3.3 0.025 0.10 0.11 130 430 15 25 1.7 430 560 140 MAIN ARM 19-24 0.6 2.5 0.033 0.088 0.10 180 400 15 25 1.7 400 580 140 25-30 0.7 2.2 0.037 0.074 0.091 190 340 15 12.5 0.8 200 390 110 31-40 0.7 1.4 0.042- 0.055 0.074 220 240 7 12.5 0.4 70 350 130 MAIN STEM 41-50 0.9 1.2 0.059 0.055 0.085 280 220 7 12.5 0.4 60 320 110 51-60 1.1 0.7 0.079 0.035 0.082 310 120 7 12.5 0.4 40 350 120 NORTH ARM 101-118 0.2 1.6 0.011 0.063 0.064 20 100 5 25 2.8 100 120 70 PITT RIVER 140-154 2.8 0.11 0.11 — 400 10 25 1.8 400 400 120 199 increase down the Main Stem - Main Arm due to increasing effects of tidal flows in the lower reaches, arid this also seems reasonable. The sensitivity of the predicted concentrations to errors in the coefficient of longitudinal dispersion is investigated in Section F.4. F.3 TIME DEPENDENT LONGITUDINAL DISPERSION COEFFICIENTS The tidally averaged values of the combined dispersion coefficient Ec listed in Table F.l represent the effects of the lateral velocity gradients on the longitudinal dispersion process when the cross-sectional mixing is essentially complete. In the initial period following the discharge of a "parcel" of effluent, the longitudinal dispersion is principally due to the effects of vertical velocity gradients alone, as discussed in Section F.l. (The time scale of vertical mixing is a half-hour). According to Fischer [1969a], the effects of vertical and lateral velocity gradients on the longitudinal dispersion are separate and additive. To allow for the variable contribution of both components, it is assumed that E = E + . E t * T 0 < t < T and } (F.10) E = E t > T c where t is the time that has elapsed since a "parce was discharged into the estuary; 1" of effluent T is the time scale of lateral mixing based on the width of section b (an edge discharge is assumed); and E , E are the coefficients of longitudinal dispersion due ^ c to the effects of vertical and lateral velocity gradients respectively. 200 Equations (F.10) are a simplistic representation of a complex three-dimen sional phenomenom, but they should adequately reproduce both the short-and long-term dispersion effects. It is only possible to account for the effects of such temporal variations in E because the tidally varying mass transport equation has been solved along its advective characteristics. Con sequently, the position of each effluent "parcel" and the time that it has spent in the estuary is known, information that is masked by a fixed grid solution to the mass transport equation. As well as increasing with time due to the influence of lateral velocity gradients as the effluent spreads over the cross-section, the coefficient of longitudinal dispersion will also vary during the tidal cycle. It will be minimal during times of slackwater and greatest during the times of strong ebb and flood flow. To investigate this effect, it was assumed that the coefficient of longitudinal dispersion given by Equations (F.9) also varied directly as the absolute velocity E = (t + S£) y U# 0 <_ t < T (F.ll) E = ay t >_ T and__ U* = °'06 5 (F.12) where E, y, and u are the instantaneous values of the respective parameters — c during the tidal cycle and a is assumed equal to the tabulated values of Ec/y in Table F.l. The relation (F.12) is obtained from Manning's equation and represents a "best over-all fit" for the entire estuary. The average of Equation (F.ll) over a tidal cycle wasrfound to be a satisfactory estimate of Equation (F.10). 201 F.4 SENSITIVITY OF PREDICTED CONCENTRATIONS The sensitivity of the predicted concentrations to assumptions about the coefficient of longitudinal dispersion was investigated for a steady discharge of a conservative effluent at Station No. 50 on the Main Stem. Figure F.2 shows the variation in concentration at Station No. 50 during a double tidal cycle and the concentration profiles along the Main Stem 50 hours after the initial discharge. The predicted concentrations have been standardized by dividing by the tidally averaged concentration ob tained from the mass of effluent discharged per tidal cycle and the fresh water discharge at Chilliwack. This is plotted as the parameter V^v/ signi fying that the results have been obtained from the tidally varying model. The effects of variation in initial dilution and multiple dosing due to flow reversal, as discussed in Section 1.3, result in spikes in the concen tration profile along the channel. It should be noted that the base of these spikes is initially only some 500 - 800 feet wide, and is much "thinner" than it appears in Figure F.2. The reason the base appears wide is that the predicted concentrations are extrapolated of the advective characteristics onto the standard 5,000 foot space grid (see Section 3.2). This grid is too coarse to accurately resolve the initial forms of the spikes. It is noted that the spikes are correctly resolved on the advective characteristics (where necessary, additional moving points were added to define regions of rapid variation, as discussed in Section 3.3). Note that the concentration of the most seaward spike (E = 0) has been halved by dilution from the Pitt River in passing through the Main Stem - Pitt River junction. Q = 36,500cfs (Chilliwack) )jan.24i|952 Tidal Range at Steveston = 11 feet J Variation in Concentration at Station No. 50 6 12 18 Time ( Hours) X tv 10 8 6 4 2 0 10 Concentration Profile along Main Arm - Main Stem (50 hours after initial discharge) (i) E = 0 (ii) E = Ey (20 ft.2/scc.) (350 ft?/sec.) t — (iii) E =E( (iv) E = Ey • Eci (v) E = (6 • oc^)y Ur 20 30 40 Stations along Main Arm - Main Stem 50 Figure F.2 Sensitivity of Tidally Varying Concentrations to the Coefficient of Longitudinal Dispersion 203 The results of assuming E to be constant and independent of time are shown as curves (i), (ii) and (iii) in Figure F.2. The corresponding values of E are: . Case (i) E = 0 Case (ii) E = E (20 square feet per second); y and Case (iii) E = Ec (350 square feet per second). The peak concentration at Station No. 50 is seen to be quite sensitive to the magnitude of E. Eighteen hours after its generation, the spike has been advected downstream to Station No. 35 and its concentration is seen to be significantly reduced irrespective of whether E equals E^ or E^. When E is allowed to vary with time according to Equations (F.10) and (F.ll), the predicted concentrations are given by curves (iv) and (v) respectively. Once again, the greatest effect is on the peak concentrations at Station No. 50, the differences between both curves being negligible after 1.8 hours when the spike is at Station No. 35. Note the significant increase in the peak predicted concentration at Station No. 50 when E is allowed to vary with u. The reason for this is apparent from the velocity variation at Station No. 50 during the tidal cycle. This is shown in Figure E.2, the spike being due to the low velocities and flow reversals around hours 4, 5 and 6. Because of the low values of u, the value of U# is very low and the effective dispersion is very small compared to other times dur ing the tidelcycle. Because of the assymetric nature of the tides, the velocities in the lower reaches of the estuary are quite low during the period of time 204 between high-high-water and low-high-water. This is illustrated by the velocity variation at Station No. 5i;in Figure E.2, and the dispersion will certainly be less during this phase of the tide than during the strong ebb and flood flow that are seen to occur once each double tide cycle. The effects of variation in initial dilution and multiple dosing are most sig nificant during the slackwaters"around the time between high-high-water and low-high-water, and this can be identified as a sensitive period of the tide cycle. In conclusion, the peak concentrations at the point of effluent discharge are very sensitive to assumptions about the form and magnitude of the coefficient of longitudinal dispersion. However, after a spike has been in the estuary several tidal cycles, its peak concentration is reason ably insensitive to the form and magnitude of the coefficient. Until adequate field data is available, the assumed temporal variations of E in Equation (F.ll) are thought to be a reasonable approximation of what occurs in the estuary. F. 5 SUMMARY The dispersion of effluent in an estuary is a complex, time-dependent three-dimensional phenomenon due to an intimate combination of the effects of turbulent diffusion, vertical and lateral velocity gradients and secondary flows. Limited field data and the results of other peoples' work have been used to obtain estimates of the coefficients of longitudinal dispersion. Due to the lack'of field data, these estimated values may be substantially in error. However, they provide a basis for obtaining preliminary notions of the tidally varying response of the estuary to waste discharges. For the assumed time-dependent behaviour of the coefficient of longitudinal dispersion the predicted concentrations were found to be relatively insensitive to the value of E^. The predicted peak concentrations at this point of efflu ent discharge were found to be very sensitive to the assumed time-dependent behaviour of the coefficient of longitudinal dispersion. It is thought that the assumed time dependent behaviour is a reasonable approximation of what happens in the estuary,-

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
Japan 7 0
France 3 0
China 2 1
United States 2 0
Canada 1 0
City Views Downloads
Tokyo 6 0
Unknown 4 13
Beijing 2 0
Ashburn 2 0
Victoria 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062998/manifest

Comment

Related Items

Admin Tools

To re-ingest this item use button below, on average re-ingesting will take 5 minutes per item.

Reingest

To clear this item from the cache, please use the button below;

Clear Item cache