Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Salinity intrusion in the Fraser River, British Columbia Hodgins, Donald Ormond 1974-01-25

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

Item Metadata


831-UBC_1974_A1 H63.pdf [ 8.62MB ]
JSON: 831-1.0064165.json
JSON-LD: 831-1.0064165-ld.json
RDF/XML (Pretty): 831-1.0064165-rdf.xml
RDF/JSON: 831-1.0064165-rdf.json
Turtle: 831-1.0064165-turtle.txt
N-Triples: 831-1.0064165-rdf-ntriples.txt
Original Record: 831-1.0064165-source.json
Full Text

Full Text

SALINITY INTRUSION IN THE FRASER RIVER, BRITISH COLUMBIA by DONALD ORHOND HODGINS B.A.Sc. University of Waterloo, 1969 M.A.Sc. University of Waterloo, 1970 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 reguired standard: THE UNIVERSITY OF BRITISH COLUMBIA AUGUST 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 Civil Engineering The University of British Columbia Vancouver 8, Canada Date 23 August 1974 ABSTRACT i The dynamics of salt water intrusion in the tidal Fraser estuary was studied by both a programme of field measurements and the use of numerical solutions of the equations of motion. Time series conductivity measurements spanning several tidal cycles indicated significant penetrations exceeding an estimated 15 kilometers above Steveston for tides of large diurnal inequality. Large ebb tides washed salt water out of the river despite low winter discharges averaging 1100 m3/sec. Mixing sufficient to disperse the salt water throughout the water column was not observed although surface currents typically ebb between 2 and 3 meters/second, and the salt wedge appeared to flood and ebb in a fairly well-defined layer. Longitudinal salinity gradients were detectable in each layer, indicating that two-way mixing took place during flood and ebb periods. Both conductivity and velocity data revealed that maximum intrusion lagged high water by 60 to 80 minutes near the river mouth. A numerical two-layer model predicted the salt water thickness within 10 per cent of the total depth and a phase agreement of ±40 minutes at maximum intrusion. Velocities were comparable to measurements within 15 cm/sec. The model neglected mixing across the interface but included the Reynold's stress formulated as KiyoD|U|Fi where Ki=0.0075, 0 is the relative layer velocity and Fi is an interfacial Froude number. The bottom stress was included as Kb^'u'lu'l where Kb=0.0055, and both stresses were found to be significant^ in the dissipation of energy in the flows. ii TABLE OF CONTENTS PAGE ABSTRACT i TABLE OF CONTENTS ii LIST OF TABLES v LIST OF FIGURES . , V NOTATION ix ACKNOWLEDGEMENTS xii CHAPTER 1 INTRODUCTION 1 The Fraser Estuary 3 Previous Work 6 CHAPTER 2 INSTRUMENTATION 10 Location of Instruments 1 Conductivity Profilers 12 Operating Principles 15 Installation and Servicing 16 Calibration 20 Data Processing and Accuracy 23 Current Metering 9 CHAPTER 3 OBSERVATIONS 34 Conductivity and Temperature 3Salinity Profiles 36 Conductivity Charts 41 Current Structure 7 Parameterization 9 iii CHAPTER 4 THEORETICAL CONSIDERATIONS 56 > The Equations of Motion 59 The Characteristic Structure and Boundary Conditions 62 A Steady Flow Model - Initial Conditions 66 The Boundary Stresses 78 CHAPTER 5 SOLVING THE TWO-LAYER MODEL ' 87 The Finite Difference Equations 88 The Barotropic Model 95 Computational Control and Boundary Conditions ...... 100 Initialization and Wash-Out 105 Comparing Predictions with Observations 106 Two-layer Results 107 Qualitative Comparisons 112 Sensitivity analysis 118 Alternate Interfacial Stress Forms 121 Essential Features of the Model 123 CHAPTER 6 SUMMARY AND CONCLUSIONS 130 Concluding Remarks 133 Recommendations for Future Research 135 BIBLIOGRAPHY 138 APPENDIX A ENGINEERING ASPECTS OF THE STUDY 141 The Effects of Mixing and Stratification .... 147 iv LIST OF TABLES TABLE PAGE 1 Comparison of derived and measured salinities for station 2 on March 29, 1973. 28 2 Summary of accuracies associated with instruments used for the Fraser River project. 33 3 Comparison of Richardson numbers obtained from Fraser River data and Taylor (1931). 82 4 Salt wedge intrusion characteristics for estuaries of various depths. 143 V LIST OF FIGURES FIGURE PAGE 1 Map of south-western British Columbia showing Fraser River and Strait of Georgia. 2 2 Fraser River hydrograph for 1973 at Agassiz, British Columbia. The Hope discharge reading is adjusted for inflows above Agassiz. 4 3 Lower Fraser estuary showing the locations of the various measuring instruments. 13 4 Cross-sectional profiles and probe line locations at the three measuring stations in the Lower Fraser estuary. 14 5 Diagram of the timber pile platforms showing general arrangement of instrumentation. 18 6 (a) Inductive conductivity probe and sensing line. Guide wire and coupling shackle appear to the left. (b) From left to right: lead acid storage batteries, Rustrak recorder, and instrument electronics. (c) Timber pile platform at station 2. 19 7 Timber pile platform at station 1 after damage. Skirt logs are missing and upper pile clamp on furthest pile is ineffective. The current speed is approximately 1.5 meters/second. 21 8 Comparison of measured conductivities by two different instruments (a,b) and measured and predicted temperatures (c,d) for two times on March 29, 1973. 24 9 (a) Operating the portable salinometer off the skiff. (b) Hydro Products current meter and trolley running along a guide wire. (c) Operating the current meter at station 2. 31 10 Typical salinity data from the Fraser estuary on a strongly ebbing tide. The Agassiz discharge is 960 m3/sec. 35 11 Conductivity and temperature data obtained with 'an RS5 salinometer on four occasions in the Fraser estuary. 37 vi 12 Salinity-depth profiles obtained on February 10, 1973 and March 18, 1973 at stations 1, 2 and 3. Channel bottom is the greatest sounding at station 1. 38 13 Contour maps of conductivity at three stations in the Fraser estuary. Bottom is marked relative to the profiler and the free surface lines are indicated using measured data at stations 1 and 3. The Agassiz discharge is 1130 m3/sec. Times are noted in Pacific Standard Time. 43 14 Contour maps of conductivity at three stations in the Fraser estuary. Bottom is marked relative to the profiler and the free surface lines are indicated using measured data at stations 1 and 3. The Agassiz discharge is 1060 m3/sec. Times are noted in Pacific Standard Time. 44 15 Current speed measurements made at two stations in the Fraser estuary on March 30, 1973. 48 16 Comparison of h1 calculated from equation (4) and the measured conductivity structure on March 18, 1973 at station 1. 51 17 Two-layer parameterization of the conductivity data obtained in February 1973. Bottom lines represent the river bed at the platform sites. Free surface lines are plotted from measured data at stations 1 and 3. 53 18 Two-layer parameterization of the conductivity data obtained in March 1973. Bottom lines represent the river bed at the platform sites. Free surface lines are plotted from measured data at stations 1 and 3. 54 19 Visualization of salinity intrusion into the Fraser estuary. 58 20 Notation used in mathematical modelling. 58 21 Possible flow states for the two-layer stratified computation. The required number of boundary conditions is indicated. 65 vii 22 Notation used in steady flow solution. 69 23 Stage-discharge relationship for Fraser estuary. 624 Measured and predicted surface elevations at three stations on the Main Arm of Fraser estuary. 74 25 Typical solutions for an arrested salt wedge. 75 26 Typical interfacial stress curves for arrested salt wedges asuming a uniform surface slope. 77 27 Profiles of °~t and velocity on March 29, 1973 at station 2. Distributions of the Brunt-Vaisala frequency and gradient Richardson number are shown. 83 28 Finite difference computation molecules. 90 29 Illustration of relationship between continuum and finite difference domains of dependence for a stable explicit difference scheme.. 93 30 Schematization of Lower Fraser River used in the barotropic hydraulic model. 96 31 Comparison of predicted surface elevations with observations at stations 3 and 4. 98 32 Comparison of predicted surface elevations with observations at stations 5 and 6. 99 33 Comparison of predicted velocities with observations on March 30, 1973 at stations 1 and 2. 101 34 Relationship between the estuary and model parameters. Conductivity profilers at stations 1, 2 and 3 are shown in their March configuration. 103 35 Comparison of typical interfacial positions predicted by the stratified flow model with two-layer parameterization of observations at stations 1, 2 and 3. 108 viii 36 Comparison of predicted and measured upper layer velocities at stations 1 and 2 for March 29 and 30, 1973. 111 37 Interfacial solutions from the stratified flow model are shown superimposed on coutour charts of constant density ( °~t) for three days in February 1973. The periods of supercritical outflow are indicated by cf <0 on the station 1 contour chart. 113 38 Interfacial solutions from the stratified flow model superimposed on coutour charts of constant density ( °~t) for three days in March 1973. The periods of supercritical outflow are indicated by c7 <0 on the station 1 contour chart. Supercritical inflow by cf >0. 114 39 Interfacial solutions for the tidal cycle on February 10, 1973. Isohalines interpolated from the data in Figure 12 are shown in each section. 116 40 Tidal variations in the interfacial and bottom stresses and the mean layer velocity shear at station 2 for March 17, 1973 (upper). Six longitudinal sections are shown for flood and ebb phases (lower). 126 41 Comparison of the measured boundary condition for h' at station 1 with the theoretical condition of Vreugdenhil (1970) (upper). Comparison of u' calculated from the stratified flow model and the theoretical relation of vreugdenhil (lower). 128 42 Advection paths of four particles released into the Main Arm of the Fraser River. The small numbers along each path indicate the time in hours following release. 145 ix Notation The following notation has been used in this thesis. Occasionally symbols have been used more than once and are subsequently listed with each meaning. a = salt water depth at x=0, a-, = coefficients obtained from a least squares analysis, A = cross-sectional area of river, A = square coefficient matrix, dimension 2n x 2n, b = storage width across river, b; = multiple regression coefficients, B = column vector of dependent variable derivatives, dimension n, c = x = dx/dt = inverse slopes of characteristic curves in x-t plane, c* = celerity of long gravity waves on a free surface, cf = celerity of internal waves on interface between two fluids of comparable but different density; e.g. water over salt water, C = constant of integration, C,Cs = conductivity in millimhos/centimeter, C = column vector of force terms, dimension n, f = column vector of force terms, dimension n, F (V) = column vector of functions of dependent variables, dimension n, FA = U/J £ gh = densimetric Froude number, F; = M/J g e 77 h*/h = interfacial Froude number, g = acceleration due to gravity = 9.800 m/sec2, h •= total water column depth in river, h« = thickness of salt water layer, i - time line in the x-t plane, j = space line in the x-t plane. X K = dimensionless stress coefficient, Ki = interfacial stress coefficient, Kb = bottom stress coefficient, Kc = inductive cell constant, L,L0 = intrusion length, m = number of characteristics entering at a given boundary, n = integer number of discrete flow layers, N (z) = (-g dp / di)/ p = Brunt-Vaisala frequency, p = mean pressure, P = solution point in x-t plane, P = solution point in the x-t plane, g = discharge per unit width (m2/sec), Q = discharge in m3/sec or m2/sec, Q = solution point in x-t plane, R = non-dimensional salt water thickness, R = solution point in x-t plane, R^ = J € gh3 / V = densimetric Reynolds number, Ri = N(z)/( du/ dz) 2 = gradient Richardson number, Ri0 = gch/u2 = overall Richardson number, Req = equivalent resistance, S = solution point in x-t plane, S = salinity in parts per thousand ), Sg = characteristic salinity of Strait of Georgia water (%«, ) , t = time, T = temperature in degrees Celcius, u = instantaneous velocity in x direction, uh = time mean barotropic velocity, xi u0 = observed layer velocity, u = time mean velocity in x direction, u* = fluctuating or turbulent velocity component in x direction, U = u-u1 = mean layer velocity shear, V = column vector of dependent variables, dimension n w = instantaneous velocity in z direction, w = time mean velocity in z direction, w* = turbulent velocity component in z direction, x,y,z = right hand Cartesian co-ordinate axes, y = Rustrak recorder ordinate, d = p~\ = specific volume, 8x = spatial increment, St = time increment, € = p' - p /p1 =Ap/p* - density contrast, 77 = h-h* = fresh water layer thichness, X =p / p* = density ratio, V = kinematic viscosity of water, P = fresh water density in gm/cm3 or Kg/m3, p% - salt water density in gm/cm3 or Kg/m3, (p-1 .0000) x10 3 = sigma-t, Ti = interfacial stress, = bottom stress, v<z = (du/dz)/-<u*w*> = kinematic eddy viscosity xii ACKNOWLEDGEMENTS The investigation into the salinity intrusion in the Fraser River has proven to be a very interesting and challenging research area and undeniably, also a frustrating one at times. Many persons have helped me by the difficult points and in particular I would like to acknowledge the guidance and valuable intuition provided by Professor M. Quick, who supervised the research. Professor T. Osborn (Oceanography) gave me encouragement during the field work and much help in converting and preparing the conductivity profilers for use. Professor P. Leblond, also in Oceanography, advised me during the theoretical studies and integrated the steady flow equations to yield equation (31). I am most grateful to both for their assistance and guidance. I also wish to express my appreciation to the Westwater Research Centre for underwriting the expenses of the field monitoring programme, and to the Institute of Oceanography for making available most of the instruments. Both organizations provided boats which were extensively used. Canada Rice Mills Limited and the Water Survey of Canada both contributed to the success of the field work and I gratefully acknowledge their help. The staff of the Civil Engineering Department were most helpful and in particular I would like to thank Mr. Geoff Sharp, whose experience in "getting things done" around the river made a valuable contribution to the field programme on more than one xiii occassion. Mr. Sharp also did a fine job in fabricating the platforms for the river. Mr. R. Brun, whose excellent work appears as the figures in this thesis is also thanked. The numerical modelling along with the data reduction were done on the IBM 370 computer available through the University of British Columbia Computing Centre and their excellent services and facilities are acknowledged. Financial support during the study period was provided by the National Research Council of Canada and the Department of Civil Engineering of the University of British Columbia. 1 Chapter 1. INTRODUCTION The Fraser River, located in south-western British • Columbia, discharges into the Strait of Georgia, a large salt water passage connected at both ends to the Pacific Ocean (Figure 1). The confluence with the Strait of Georgia results in a salt water intrusion along the river bottom extending several kilometers above the mouth. Observations have shown that in many respects a circulation pattern typical of salt wedge estuaries is established in which salt water flows landward along the bottom while river water flows seaward in a relatively fresh surface layer. Unlike other coastal plain estuaries where well defined salt wedges can be maintained, the very large tides produce intensified flows of tidal period in each layer. The results of a study into these tidal circulations carried out during 1972 and 73, are reported in this thesis. The main objective of the project was to examine the nature of the seawater intrusion and arrive at a quantitative relation between the Strait of Georgia tides and the observed currents and density structure. This meant that several related problems had to be tackled along the way, such as arriving at a description for the stress along an interface between two fluids of different density. The study proceeded along two main directions. First, a programme of field measurements was undertaken during the winter of 1973 to collect suitable data to provide a detailed picture of the salt wedge movements on a sub-tidal scale. To complement, 3 and help explain the field observations, mathematical models of stratified tidal flows were developed. The Fraser Estuary The lower Fraser River, which we may define as extending from Hope to the Strait of Georgia (Figure 1), undergoes relatively little change in elevation resulting in flows that are tidal over as much as 110 kilometers during low discharges. Since the river has an extensive watershed and a significant proportion of the flow is derived from snowmelt, there is a large seasonal flow variation. The annual discharge curve for 1973, shown in Figure 2, is typical and the low winter flows result in the density stratified estuary conditions which are the subject of this study. The lower estuary divides into two principal channels and a number of lesser branches, some of which are connected to the Strait of Georgia only at high tide. This investigation has been confined to the main arm of the river which carries approximately 80 percent of the,flow. It is also the deepest channel with an average depth of 10 to 15 meters throughout the seaward 30 kilometers. The main arm is continuously dredged to provide a minimum shipping channel depth of 9 meters, a feature which has probably accentuated the salinity intrusion. An understanding of the tidal hydraulics in the estuary is of fundamental importance. There are, perhaps, two main concerns: sedimentation and water quality. Preserving the 4 Figure 2. Fraser River hydrograph for 1973 at Agassiz, British Columbia. The Hope discharge reading is adjusted for inflows above Agassiz. 5 shipping channel is a year round dredging operation, consuming large public expenditures, and made necessary by continuous sedimentation of the water-borne silt load. It is commonly observed that deposition is accentuated in the vicinity of salt wedges; therefore, a predictive capacity for the stratified river hydraulics would be of value in the sedimentation studies. The lower estuary is flanked on one side by the city of Vancouver and a number of smaller communities on the other. As a result, the concentrated use of the river, principally for transportation, timber movement and storage and domestic and industrial waste disposal, has recently produced concern about the water's quality. Assessing water quality, either by using numerical dispersion models or field monitoring is complicated by the salt water intrusion and relies on a knowledge of the principal circulation patterns. The field investigation forming part of this project, was the first intensive effort to gather data on the tidal variations of the salt wedge motion in the river and a description of the instruments, some of the techniques and problems is given in Chapter 2., Some typical examples of the data along with a brief discussion of the significant features is presented in Chapter 3. In Chapter 4 the basic mathematical models are derived, and an analytic solution for a steady flow model is presented. The numerical procedures are developed and applied in Chapter 5 together with a comparison of observed and predicted results. The study is ^concluded in Chapter 6 with a discussion of some 6 important aspects of salinity intrusion in the estuary. Previous Work Considerable research effort has been directed toward controlling sedimentation in the lower Fraser over the last twenty years. Notable in this regard was the physical scale model constructed at the University of British Columbia by Professor E. Pretious on which various channel controlling schemes were developed. Two extensive series of stage measurements and the cubature calculations of discharge which can be made from these data, are discussed by H. D. Baines (1952,1953). The series were made for high tide-low flow and low tide-high flow conditions in 1952. However, research during this period was confined to unstratified river conditions and no attempts were made to relate the hydraulics and sedimentation to the salt penetration. Although marine borers in timber piling have long indicated salt water presence in the estuary, physical measurements of the density structure for various tidal conditions were rarely made. One series of salinity/temperature/depth measurements was reported by Waldichuk, Markert and Meickle (1967). Several series of salinity and temperature profiles at various locations in the delta were made during 1971 by A. Ages of the Marine Sciences Branch of Environment Canada (unpublished) and these data were valuable in illustrating stratified conditions for high tides 7 during the winter months, and aided in the design of the instrumentation. The Water Survey Branch of Environment Canada continuously monitors the freshwater discharge in the Fraser River at Hope, B.C. (Figure 1), and the river surface elevation at eight locations upstream of the Strait of Georgia. The Department of Public Works also publishes very detailed semi-annual sounding surveys made between the Strait of Georgia and Pitt River. These are the only physical data for the Fraser estuary available on a continuous basis. The hydraulic mechanisms controlling salt penetration in "tideless" estuaries have been studied theoretically by Keulegan (1966) and Farmer and Morgan (1953). Keulegan applied his results for arrested salt wedges to the Mississippi River where success in predicting penetration lengths was achieved. Farmer and Morgan obtained good comparisons between a "steady-state" model and experimental flume data. However, it is clear that the Fraser estuary is far from a steady flow and the tidally induced motions are important. From the standpoint of mathematical modelling, time varying equations are required and numerical solution offered a good chance of success. Two numerical models describing the unsteady behaviour of salt wedges have been published by Vreugdenhil (1970) and Boulot, Braconnot and Marvaud (1967). In both analyses a two layer model was assumed and conservation of mass and momentum equations for each layer were 8 numerically integrated using finite difference schemes. Boulot et al. did not include salt water exchange between the layers in their model and it can be classified as a no-mixing analysis. They presented results showing the salt front oscillating over a distance of less than one kilometer. This short distance resulted from the small tides which averaged 0.14 meters in amplitude. The model was applied to the Grand-Rhone downstream of Aries, but unfortunately no correlations between their model predictions and physical measurements were discussed. Vreugdenhil*s two layer model also neglects mixing of salt and water between the layers, and results from the model are discussed in relation to observed prototype measurements from the Rotterdam Waterway in The Netherlands. Using essentially the same equations as Boulot et al., Vreugdenhil's method succeeds in reproducing the main salt water intrusion features. Of interest is the salt wedge shape shown for several times in the tidal cycle. The mean interfacial slope ranges from 3x10~4 (ebb) to 8x10-4 (flood) (estimated from data shown) suggesting a considerable variation in the salt water velocity. Penetration distances vary between 9 and 2 Km for the 1.6 meter tide. Vreugdenhil concludes that the empirical coefficient which enters the analysis in the interfacial stress term has a dominant effect on the model and states he was able to calibrate the predictions to give a good correspondence with observations using this coefficient. By using a calculation for the horizontal tide, that is 9 the barotropic or unstratified velocity field as a function of time and position in the estuary, Vreugdenhil reduced the total number of equations required for the simulation from four to two. As a result he found that the overall behaviour of the salt wedge was largely dependent on the horizontal tide prediction. 10 Chapter 2. INSTRUMENTATION Previous field measurements from the Fraser River -indicate the presence of salt water along the river bottom during winter months. However, variations in the extent of salt water penetration and correlations between the salinity structure and Strait of Georgia tides are difficult to determine from these data. Tidal amplitudes at the river entrance can exceed 40 percent of the depth and are sufficient to reverse the river flow for several kilometers above the region of seawater penetration. Under these conditions, the value of point measurements in time is limited. It was clear that time series data would be reguired to adequately describe the density field and a programme of measurements in the river, emphasizing continuous records of the conductivity structure, was undertaken during the winter of 1972/73. Three instruments described by Farmer and Osborn (1973) were installed during February and March of 1973 which provided conductivity over most of the water column every 15 minutes,"and are subsequently described as "conductivity profilers". The actual periods of installation were: February 01 to 13 and March 16 to 30. The two most important forcing functions for the Fraser River are the tide and the freshwater discharge. Both of these parameters are continuously monitored by the Water Survey of Canada (Environment Canada). The discharge is computed from a calibrated gauging staff at Hope, B.C. on an average daily basis. 11 This figure is then adjusted for Agassiz, B.C., the furthest upstream point of tidal influence. The Agassiz figure will be used throughout. During the final two days of March, current metering was carried out simultaneously with conductivity measurements. This provided comparable velocity and salinity data. At other times during the winter, salinity/temperature/depth (STD) measurements were made using an Industrial Instruments RS5 salinometer. The sensing head was attached to a winch line above a 25 pound weight and operated off the stern of a small boat. The winch used was capable of very accurate depth readings but the real limitation to accurate depths was the current speed. Ebb currents produced large wire angles (exceeding 20 °) with this light equipment over a considerable portion of the ebb phase and measurements were abandoned for angles exceeding 30 °. This illustrates the main difficulty in making physical measurements in the Fraser River, and helps to explain why most of the available data were obtained during high water slack periods - occasions when salt water was well established in the river. Location of Instruments In any scheme to collect time series data from estuaries many factors influence instrument disposition. Securing and protecting the conductivity profilers in the large ebb currents, while not interfering with navigation facilities, were the main considerations in the Fraser River. .It was also important to 12 monitor as much of the water column as possible, and locate the profilers adjacent to the government maintained tide gauges. Figure 3 shows the Fraser River with the enlarged chart, covering part of the region of salt water penetration, indicating the locations of the various permanent instruments. Conductivity profilers 1 and 3 were located adjacent to the tide gauges, and profiler 2 was placed between these two. Timber pile platforms were erected in the river beside the navigation channel for profilers 1 and 2. This provided a rigid system for mooring the sensing lines and a degree of protection from debris, river traffic and currents. At low tides, a minimum thickness of 6.1 meters of water was monitored. The third profiler was located on the pier of Canada Rice Mills Limited, where two separate installations were made to obtain improved coverage of the water column. The cross-sectional shape of the river and the extent of coverage for each profiler, are shown in Figure 4 for stations 1, 2 and 3. Conductivity Profilers The salinity and temperature data collected during the winter months from the Fraser River illustrate some characteristic properties of the two water masses. During those tidal periods when salt water from the Strait of Georgia was established along the river bottom, salinities typically ranged from 25 to 29 parts per thousand, and temperatures of 6 to 8° C were found. The fresh water had temperatures ranging from 2 to Figure 3. Lower Fraser estuary showing the locations of the various measuring instruments. 14 G.D. LOW WATER LEVEL J STATION 12 G.D. LOW WATER LEVEL STATION 2 \ G. D. |0 2 4 6 8 10 METERS STATION 3 100 200 METERS Figure 4. Cross-sectional profiles and probe line locations at the three measuring stations in the Lower Fraser estuary. 15 5° C throughout the measuring period and no detectable salinity. For these temperature ranges both the conductivity and density profiles are primarily a function of the salinity structure. Thus a measure of the conductivity distribution can be used to infer the salinity profile and hence the density structure. Ideally, both conductivity and temperature should be recorded; however, these estimates can be improved by incorporating the temperature data available from the STD measurements. Operating Principles The instruments measure the electrical conductivity of sea water at fourteen depths by successively interrogating inductive probes on a sensing line secured in the water column. A fifteenth probe, located inside the instrument housing, is also interrogated in every cycle to act as a check on any drift in the electrical response. Each probe consists of two mu-metal annular cores surrounded by torroidal coils of insulated copper wire. A constant voltage is applied to one coil, which induces an electrical current in its core. This current in turn induces an electrical current in the sea water path linking the two coils. The second coil acts in a complementary fashion to the first, carrying a current proportional to the electrical conductivity of the sea water path. Each probe was interrogated for one minute and the output recorded on a Rustrak chart recorder. The sensing line was 16 interrogated from bottom to top, with the calibration probe occupying the fifteenth position. In this way, a sequence was completed every fifteen minutes, producing four cycles to the hour, a frequency considered satisfactory for examining tidal phenomena. Timing of the interrogation interval was controlled by a small hobbyist's clock, operating on an independent power supply inside the instrument. These clocks, with an accuracy exceeding ±4 minutes per week, made the timing independent of the Rustrak behaviour, a desirable feature for outdoor operating conditions. The instruments were powered by four 6 volt lead acid storage batteries, which provided a near constant supply voltage for two weeks. Rustrak charts required replacement once each week. Installation and Servicing Ebb tides regularly produce currents exceeding 2 meters/second in the Fraser River. In addition, the river carries considerable boat traffic of all types, and floating debris primarily from the logging industry. Holding and protecting a sensing line securely in place, under the constraint of reasonable cost, provided a real challenge. A platform seemed to be the best solution because it could serve a variety of purposes. Servicing from a small boat would be possible, and other types of measuring equipment could be operated in conjunction with the conductivity profilers. Any platform also had to be free of strong oscillations which would 17 affect the Rustrak operation and could lead to structural deterioration. A final design of three timber piles arranged in an equilateral triangle was chosen. Two steel tube frames spaced three feet apart ^joined the piles. These frames were crossbraced and the lower triangle supported a working deck. The sketch in Figure 5 gives the pertinent dimensions and arrangement of the instrument and Figure 6 (c) shows the actual installation at station 2. The perimeter triangle of logs, floating just outside the piles, provided some protection for the probe line by deflecting debris. The sensing line was lowered down a guide wire secured to the platform, stressed by a heavy weight resting on the bottom (Figure 5). The probe line was coupled in five places to the guide wire by brass snap-shackles, often used in sail boat rigging. The inductive probe and sensing line, coupled to the guide wire at the left, are shown in Figure 6 (a). A 30 pound weight on the probe line kept it extended. s Once installed, the instruments were easy to service. The platform was large enough for one person to raise and lower the probe lines conveniently, and batteries could be exchanged without difficulty. As shown in Figure 6 (b), the batteries, instrument and Rustrak were arranged in a weather-tight plywood box. This meant that Rustrak charts could be replaced under cover during bad weather, a frequent occurrence in winter. The plywood box also provided some security against vandalism. 18 GUIDE WIREN PROBE LINE. FLOOD CURRENT EBB CURRENT C C C C C O PLAN VIEW HIGH WATER LOW WATER PROBES (14) GUIDE WIRE -SLIDING WEIGHT -PERMANENT WEIGHT ELEVATION VIEW Figure 5. Diagram of the timber pile platforms showing general arrangement of instrumentation. 19 Fijure 6. (a) Inductive conductivity probe and sensing line. Guide wire and coupling shackle appear to the left, (b) From left to right: lead acid storage batteries, Rustrak recorder and instrument electronics. (c) Timber pile platform at station 2. 20 Two difficulties were encountered. The ebb currents were so strong that significant probe line deflections occurred during the first two week period of installation. To alleviate this problem more weights were added to the guide wires, increasing the total submerged weight to 185 pounds. For future installations of this type a minimum figure of 300 pounds is recommended. During the last week of February platform 1 (station 1) was damaged. The skirt logs were smashed free of the piles, causing the entire structure to go approximately 5 degrees off vertical. A weld on the upper triangle pile clamp was also broken and the resulting loss in rigidity made the platform more susceptible to flow-induced oscillations. Consequently, progressive deterioration of the remaining welds occurred over the next two weeks, until one third of the securing clamps were ineffective. Despite this and the increased amplitude of oscillation, the platform remained operative to the end, attesting to the design's effectiveness. Figure 7 shows platform 1 after sustaining damage. The current is ebbing at approximately 1.5 meters/second in the photograph. Calibration Following the procedure developed by Farmer and Osborn (1973), the calibration of each probe was carried out in two stages. First a cell constant Kc, depending only on the geometry of the probe casing, was found. The value of this constant 21 Figure 7. Timber pile platform at station 1 after daiaage. Skirt logs are missing and upper pile clamp on furthest pile is ineffective. The current speed is approximately 1.5 meters/second. 22 varies slightly from probe to probe but it is always independent of the solution conductivity and the electrical response of the instrument and recorder. The final stage involved calibrating the electronics and chart recorder as a function of conductivity. The cell constants were determined by immersing each probe in several sea water solutions of known conductivity. The output voltage was recorded for each solution, after which the probes were removed from the water and dried. A conducting loop in series with a variable resistance box was passed through the opening in each probe casing, and the "equivalent resistance", Reg, found which yielded the same output voltage obtained in the sea water solution. Knowing the solution conductivity, Cs, the cell constant is defined by: Kc = Reg. Cs (1) Three separate solutions of increasing conductivity were used to find Kc and the average value of these three was used for each probe. The worst scatter in the averaging process was typically about 11 percent of the mean value. Once the cell constants were known, the chart recorder and electronics were calibrated by simply adjusting the resistance box through a sequence of resistances compatible with the conductivities to be measured in the field and noting the recorder output. In practice, nine resistance values were selected ranging from 100 to 1000 ohms. A least squares analysis 23 was used to find a functional approximation between the recorder output and the resistances. The final equation for reducing the chart record has the form: conductivity (mmhos/cm) = Kc a, +a2y + a3y: (2) where a-, are the least squares coefficients and y is the recorder output. Data Processing and Accuracy, Farmer and Osborn (1973) quote an absolute accuracy of ±2 mmhos/cm and a relative accuracy of ±1 mmho/cm between the probes on any one line of the conductivity profilers. These values apply to a range of conductivities of 0 to 35 mmhos/cm. I found the repeatability of individual measurements during calibration to lie within the +2 mmho/cm figure and have applied this value to the results discussed in subsequent sections. In Figures 8 (a) and 8(b) nearly simultaneous profiles of conductivity obtained from two different instruments are plotted for two times at j i station 2 on March 29, 1973. The shaded band on each conductivity profiler line represents the range of possible values based on an accuracy of ±2 mmhos/cm. The values obtained from the RS5 salinometer are shown with an error bar of ±1 mmho/cm. The profiles were recorded within 15 minutes of one another for each comparison and the results in Figure 8 (a) illustrate the worst field comparison obtained, and those in 8(b) one of the best. In the'latter case the accuracy of ±2 mmhos/cm 24 CONDUCTIVITY ( MMHOS/CM ) TEMPERATURE (°C ) STATION 2 MARCH 29,1973 • CONDUCTIVITY PROFILER KH RS5 SALINOMETER Figure 8. Comparison of measured conductivities by two- different instruments (a,b) and measured and predicted temperatures (c,d) for two times on March 29, 1973. 25 seems very reasonable, as it does for both the near surface and two deep probes in Figure 8(a). The profiler results cn either side of the line in 8 (a) are almost identical in shape and magnitude, suggesting a structure that is quite stable following a small ebb tide. Therefore, the departure of the profiler results from those of the salinometer between 2 and 5 meters is due either to a shift in the water column structure between measurements or calibration factors that are in error for certain j probes. The second cause is most likely and a similar trend can be observed for the same probes in Figure 8 (b) although the. differences are smaller. Since the ±2 mmho/cm accuracy on the profiler results always overlapped the RS5 values (including their error), this accuracy figure was considered applicable. Two factors seemed to influence the calibration values for each probe: humidity and temperature effects on the instrument electronics and moisture seepage and condensation inside the probe heads. The exact effect of humidity and temperature changes on the electronics was not determined but in order to minimize their influence calibration was done under outdoor conditions. . A more serious problem concerned moisture effects within the probe heads. Often the deep probes would become flooded, fail entirely and were easily detected. In other cases a few beads of moisture were apparent after removal to the laboratory but the appearance of the Rustrak charts was normal. In practice these "moist" heads, as distinguishable from flooded probes, were found for three probes during the final installation period and their data were retained in the analyses. The data 26 presented in subsequent sections were reduced from the Rustrak charts using calibration factors obtained immediately prior to installation and no attempt was made to account for moisture induced changes. Probe heads were eliminated if failure could be positively determined. Salinity and density data were derived from the measured conductivities and used in the parameterization of the results. A reasonable estimate of salinity can be made without including temperature effects but incorporating background temperature data improved the accuracy of the estimates. Temperatures were computed for each conductivity based on linear relations between the conductivity and temperature data accumulated from the RS5 salinometer surveys for particular installation periods. Comparisons of measured and derived temperatures are shown in Figures 8(c) and (d) , where an accuracy of ±0.5 °C would be representative of data derived in this way. Knowing conductivity and temperature, the salinity was computed from the following relation: S = b, C + b2T + bjCT, (3) where S = salinity (parts per thousand), C = conductivity (mmhos/cm), T = temperature (°C), as derived from the conductivity data. The coefficients bj were obtained from a multiple regression 27 analysis performed on values of the three parameters obtained from Lafond»s Tables (Lafond (1951)). Using the field profiles from Figures 8(b) and (d), computed and measured salinities are compared in Table 1. An accuracy is listed for each computed salinity based on equation (3) and tho previously quoted error ranges for conductivity and temperature. The differences between calculated and measured quantities appear in column 6 and have a root mean square value of 0.53 parts per thousand. It appears, then, that the salinities derived from the conductivity profiles are considerably closer to the actual values than the computed error averaging ±3. 1 %o would indicate. In fact the error range of ±3.1 %o represents the maximum limit if all inaccuracies in the conductivity and temperature data combined in the worst possible way. Density was computed using the relation given by Lafond (1951) in terms of sigma-t, °~t, which can be calculated to an accuracy of ± 0.01 for given values of salinity and temperature. Sigma-t is related to the specific gravity of sea water, p , by °~t = (P -1.0000) x103 and properly diroensionless. Since the density of fresh water is taken as 1.0000 gm/cm3, the specific gravity of sea water is often written in terms of density (that is given units of gm/cm3 or Kg/m3) and this has been done here. To estimate an accuracy for O" t calculated from the profiler results, conductivity and temperature data including their errors were combined in the worst way and evaluated in Lafond's equation. The maximum range was established and one-half of this figure used for the accuracy. Column (7) of Table 1 lists a Table 1 Comparison of derived and measured salinities for station 2 on March 29-, 1973 (1552 hours). * indicates interpolated values. SALINITY DEPTH COND. TEMP. CALCULATED MEASURED* DIFFERENCE SIGMA-T (meters) (mmhos/cm) (°C) (7oc) (%o) (%o) (D (2) (3) CO (5) (4)- (5) = (6) (7) 0.3 4.47±2.0 5.89±0.5 3.83±2.9 3.00±0.5 + 0.83 3.04±2.3 0.8 3.89 5.86 3.24±2.9 3.25 -0.01 2.58±2.3 1.3 4.61 5.89 3.97±3.0 3.60 + 0.37 3.16±2.3 1.8 5.01 5.91 4.37±3.0 4.30 + 0.07 3.47±2.3 2.8 6.65 6.00 6.02+3.0 6.75 -0.73 4.77±2.3 3.3 9.06 6.12 8.42+3.1 9.00 -0.58 6.66±2.4 3.8 12.09 6.27 11.43±3.1 11.20 + 0.23 9.02±2.4 4.3 14. 18 6. 38 13.49±3.2 14.15 -0.66 10.63±2.5 7.3 . 26.34 6.99 25.28±3.4 - - 19.8U2.7 RMS value 0.53 NJ CO 29 sequence of °"t • s calculated from the profile in Figure 8(b) together with the errors. An average value of ± 2.4 was obtained and as in the case of salinity this figure represents the maximum range. C«££gnt Metering One aspect of the Fraser River seldom examined by field measurements is the current structure and its relation to the tide. Due to large wire angles, making reliable current measurements from a small boat is virtually impossible once the river velocities exceed one meter/second. In addition, the boat motion itself can induce large errors, particularly in the current direction, if the meter is suspended from the vessel. However, once the platforms were installed it became feasible to hold a current meter in the flow at a constant wire angle. The Water Survey of Canada has operated a catamaran with current speed measuring equipment on the Fraser River for several years. A joint velocity measuring programme was undertaken on March 30, 1973 with the Water Survey boat occupying station 1 and our equipment located at station 2. The object was to obtain a siumlataneous record of velocity and density as a function of depth and tide. At this time both conductivity profilers were functioning at these stations; however, three probes at station 1 and two at station 2 were lost. The platform at station 2 was used as a fixed point to suspend the current meter. To eliminate the problem of changing 30 wire angles with changes in current, a guide wire running from the platform to a 300 pound anchor was installed. This gave a constant wire angle of 10 degrees in all currents. The current meter was equipped with a trolley running on the guide wire (Figure 9(b)) and this system was raised and lowered by hand using a calibrated rope to give the required depth. To obtain the most stable configuration, the hauling rope was tied off to the platform for each reading. Figure 9 (c) shows the current meter and the installation and in Figure 9 (a) the method of operating the RS5 salinometer is shown. The current meter was a Hydro Products 400 series speed and direction meter. A Savonius rotor activates the speed sensor, and a solid vane provides direction readings relative to a magnetic compass. An accuracy of ± 3 percent for the speed and ± 5 degrees in the direction is quoted by the manufacturer. In general, the direction values were stable for each depth except when the speeds became small. The guide wire had a minimum clearance of 1 1/2 meters from the nearest pile and was located in a manner to minimize platform interference of the flow. The Water Survey of Canada eguipment was designed primarily for estimating discharges in unstratified rivers. The current meter is suspended on a single line just above a 150 pound weight. The magnitude of the current vector is obtained but the direction relative to the river is not measured. The channel at station 1 is fairly straight and the principal flow 31 Figure 9. (a) Operating the portable salinometer ofr the skiff. (b) Hydro Products current meter and trolley running along a guide wire. (c) Operating the current meter at station 2. 32 directions are more likely to be parallel to the shoreline than at station 2. Although the actual velocity component in a desired direction cannot be computed from the station 1 data, it is reasonable to assume that current speeds exceeding 20 cm/sec represent flows aligned with the banks. Station 2 was also occupied during daylight hours on March 29, 1973. Velocity profiles were measured every 45 minutes, as on March 30. A summary of the accuracies associated with each instrument and the derived data for the Fraser River project is presented in Table 2. 33 Table 2 Summary of accuracies associated with instruments used for the Fraser River project. * indicates derived values as explained in the text. CONDUCTIVITY RS5 HYDRO PRODUCTS PROFILER SALINOMETER CURRENT METER (D (2) (3) CONDUCTIVITY (mmhos/cm) ±2.0 ±0.5 -TEMPERATURE (°C) ±0.5* ±0.5 -SALINITY < %o> ±3. 1* ±0.3 -°t ±2.a* - -SPEED (cm/sec) — - ± 3 PERCENT DIRECTION (° of arc) - - ±5.0 34 Chapter 3. OBSERVATIONS On five occasions during the study period salinity/temperature/depth surveys were made. These data were collected primarily to act as a field check on the conductivity profilers and to provide background information on the temperature structure. Making these profiles was fairly quick once the boat was on station and allowed an examination of the salt water structure at various locations in the river over a short space of time. On March 8, 1973 four stations were sampled in the space of 60 minutes, spanning a distance of approximately 10.5 kilometers. The first two stations coincide with conductivity profilers 1 and 2 and the upstream locations straddle station 3. At the fourth station, located in the deepest part of the channel, no traces of salt water were found. These data are profiled in Figure 10 and a longitudinal section is plotted by' interpolating for the isohaline depths at each station. One obvious feature is the wedge like shape of the salt water intrusion. Despite the strong ebb currents - profiling at buoy 28 was done at limiting wire angles - the stratification is fairly well maintained, especially between stations 1 and 2. It can also be seen that the average salinity in the salt wedge decreases in the upstream direction indicating that fresh water mixes downward to dilute the salt water below. Conductivity and Temperature The measured conductivities reflect the extent of mixing 35 SALINITY (%o) SALINITY (%o) (d) 2 BUOY 27 BUOY 28 0 2 4 6 8 DISTANCE (KM) (e) Figure 10. Typical salinity data froin the Fraser estuary on a strongly ebbing tide. The Agassiz discharge is 160 m3/sec. 36 of salt between the two water masses and as a rough approximation, the temperatures can be assumed to be distributed in a like manner and linearly related to conductivity. The temperature data collected during STD measurements are summarized in Figure 11 and if the data is considered in three groups separated from one another by roughly 30 day intervals, a linear relation between temperature and conductivity for each group is very reasonable. Since both water masses warmed during the measuring period, the pertinent temperature relation from Figure 11 was used to calculate salinity and density from the conductivity profiler results for a particular two week installation as described in the last section of Chapter 2. This type of calculation gives a reasonable estimate considering the overall accuracy of the conductivity profilers and the effect of temperatures in this range on salinity and density. Salinity Profiles There is a considerable variation in the diurnal inequality of the Fraser estuary tidal amplitudes, ranging from nearly zero when each succeeding tide looks much as its predecessor to very large, which produces a high-low water only a fraction of a meter below high-high water. This variation produces markedly different intrusion characteristics as revealed in the conductivity data. In Figure 12 two series of salinity profiles have been plotted which illustrate the behaviour at stations 1, 2 and 3 for each extreme tide type. The salinities were derived from the conductivity profiler data using equation 37 8 STATION I A JAN. 30 0925 O JAN. 30 I 300 O MAR. 8 1039 STATION 2 © JAN. 29 0920 + JAN. 29 1250 • MAR. 8 I I 0 I 0 MAR. 29 0853 9 MAR. 29 I 552 ~0" o BOTTOM WATER SURFACE WATER ± 0 10 15 CONDUCTIVITY ( 20 MMHOS/CM) Figure 11 obtained occassions Conductivity and temperature with an RS5 salinometer on in the Fraser estuary. 25 da ta four 30 Station I River discharge M30mVsec. a> 4 E 20 22 0 2 4 6 B 10 12 14 * 16 18 20 22 0 Hour* February 10, 1973 Salinity [%o 10 20 >{k 0355 H°"'» l' 2 • Station 1 • Station 2 * Station 3 •0 20 10 10 - 20 K IQ 2Q Salinity (%o) •/ 6 VI 0856 ^ \ 0933 \ V \ \ \\ \ • Station 1 * StOtion 2 * StoNttn 3 Station I River discharge l060mVsec. •'Deepest Channel bottom at Station I 22 0 2 6 8 10 12 14 16 March 18, 1973 20 22 0 rtoun Salinity (%„) *. u *\ 0258Hou't * > tt ' " ff 0358 I \ X \ \f 0456 } \ f .6 •U 10 ZD r : 0658 "'V V\ HJ 10 20 ; Y\ * 07 SB \ ! \ •0 10 20 5 \ \ • " i V i 6 io 20 3 \ "V ' ^ 09S9 { B O 20 3 1 \ ' * \ \ > 1 i i \ » j|* • S'01'On 1 t* * Siot.on 2 \ \ i \ V \ \ | • SiO'.on J Deepest Channel bottom ot Station I Figure 12. Salinity-depth profiles obtained on February 10, 1973 and March 18, 1973 at stations 1, 2 and 3. Channel bottom is the greatest sounding at station 1. CO CO 39 (3) and the linear temperature relations of Figure 11. Each profile is plotted with respect to the free surface at the particular station, not a horizontal datum, and I estimate this could introduce a maximum error of approximately 50 centimeters between the profiles for stations 1 and 3. The error would be much less for most of each series since the measurements occur nearer to high water. The solid line below each profile represents the deepest channel bottom at station 1 and as can be seen, the combination of lost probes and platform position has meant that only 70 percent of the water column is monitored at any station for the first series. The situation is slightly worse for the March series where an additional probe has failed at station 1. This means that we are watching the behaviour of the fresh water layer and seeing little of the actual salt wedge motion until the intrusion has reached at least five or six kilometers above station 2. Furthermore, station 1 is the shallowest with respect to datum of any station so this problem is worse at the upstream positions. Despite these limitations there is much to be gained by s examining the data in this manner. In the February series the long high water duration has resulted in a very long salt wedge penetration, probably exceeding 15 kilometers above station 1. The profiler at station 3 has clearly measured the distinct halocline associated with the salt-fresh interface around high tide and at all three stations salinities exceeding 25 %•> have been recorded at maximum penetration. This pattern is contrasted by the March series, where the shorter high water duration has 40 produced a significantly shorter penetration and profiler 3 seems to be measuring near the intrusion limit. In spite of the improved coverage (profiler 3 was relocated between series to gain nearly 3 more meters of water column) a less distinct halocline is evident and two-way mixing of salt and fresh water is suggested by the profiles following high tide (17 and 18). The profiles recorded at station 1 for both series show that just after high water, surface salinities reach a maximum, typically between 10 to 15 %o. In the February data any distinct halocline has been lost to an almost uniform increase of salinity with depth at stations 1 and 2. The halocline is preserved in the March series and profiles 14 through 17 for instruments 1 and 2 suggest a well mixed surface layer flowing over the salt water layer. These high surface salinities could result from strong local mixing between the water masses or the convection of salty surface waters back upstream by the flooding tide, or a combination of these mechanisms. The flood tide duration is sufficient to convect brackish surface waters, formed by mixing during the previous ebb tide, upstream as far as stations 1 and 2 from the vicinity . of Sandheads (near the final confluence of river and Strait of Georgia ). These profiles certainly show that the Fraser estuary does not have the high localized salinity gradients associated with arrested salt wedges - tidal mixing is strong enough to produce a more dispersed density structure and there is a measurable longitudinal salinity gradient in the surface and 41 bottom waters near high tide. However, the behaviour of the estuary can be viewed in terms of two layers which are discernable in these data throughout most of the tidal phase. The changes in motion in each layer will be a response to the particular force balances in the layer, comprised of pressure and friction. If the exchange of mass and momentum across the salt-fresh interface is large then the layer motions will be highly coupled and the relative phase separation of the flows will be small. Conversely if the exchange is small the relative motions will have a large separation in time. In the February data the maximum penetration of salt water occurs about one to two hours following high tide at station 3, and just after two hours at stations 1 and 2. Profiles 9 and 10 show that the vertical salinity (and density) gradients are increased during a period when the surface water is being accelerated downstream due to the seaward pressure forces. This suggests that the salt water at depth continues to flow upstream. The March series confirms these observations, where the maximum penetration seems to occur about three hours following high tide. In the absence of velocity data during these measurements it is difficult to quantify the phase relationships but it appears that the maximum penetration of salt water lags high tide at station 1 by approximately two to three hours. Conductivity Charts The conductivity profilers were operated basically for two week periods, since by the end of 13 or 14 days in the river 42 sufficient damage occurred to the probe heads at stations 1 and 2 to warrant removal and repair. A convenient method of presenting these time series data is in the form of a contour map of constant values of conductivity obtained by interpolating between the measured values for each cycle. Conductivity and tidal data for two different types of tide are plotted in Figures 13 and 14 where the vertical axes for the contour maps represent the distance in meters above probe number one, and the horizontal axis is time. The bottom and free surface lines are also plotted. Tidal heights above Sandheads datum were used directly for stations 1 and 3 since the probe line position relative to this datum is known. At station 2 a broken line indicates the surface based on tidal data from station 1. The depth to bottom from probe 1 was estimated from the anchor dimensions and unlike the bottom lines in Figure 12, these represent the river bed at the platform locations. As explained previously probes would fail at any time during the installation period. The time of failure could only be determined by examining the Rustrak charts and reduced data. One assumption was applied to the data: in general the water column would be statically stable within the error bounds of the individual measurements. To indicate which probes were operative for a section of contours zeros were plotted, at the appropriate depths, to the left of the section. Records . of surface height at station 1 relative to the geodetic level datura (G.D.) are also plotted in Figures 13 and '1 STN.I FEBRUARY 1973 Figure 13. Contour maps of conductivity at three stations, in the Fraser estuary. Bottom is marked relative to the profiler and the free surface OJ lines are indicated using measured data at Stations 1 and 3. The Agassiz discharge is 1130 m3/sec. Times are noted in Pacific Standard Time. STN.I MARCH 1973 Figure 14. Contour maps of conductivity at. three stations in the Fraser estuary. Bottom is marked relative to the profiler and the free surface lines are indicated using measured data at stations 1 and 3. The Agassiz discharge is 1060 mVsec. Times are noted in Pacific Standard Time. 45 14. Chart datum for the Fraser River is a sloping line representing the lowest tidal height and is therefore a function of position. Geodetic datum is considered a level plane and as such a more convenient reference datum. The contour charts in Figure 13 were obtained from the data recorded between February 8 and 11, 1973 for tides having a pronounced diurnal inequality. In contrast the tides which occurred from March 16 to 20, 1973 have a small diurnal inequality and produced the conductivity data contoured in Figure 14. Care must be exercised when interpreting the February data. As mentioned previously, probe line deflections giving a wire angle of approximately 20 degrees occurred at stations 1 and 2 during the first installation period on the large daily ebb. Observations during the remaining tidal phases indicate that excessive wire angles did not occur. Probe one at both stations 1 and 2 did not function for the entire recording period, and taking this into account I estimate that during the large ebb the lowest points a of measurement correspond to a reading of approximately 2.5 to 3.0 meters on the scale shown. The effect of this error on the contours (no attempt has been made to adjust the data) is to indicate an absence of salt water for longer periods than is the case. This problem was corrected for the March data where wire angles were confined to less than 10 degrees. The conductivity profiler at station 3 produced data for 46 only one tidal cycle in Figure 13 (seven probes), but was operative for the entire period in Figure 14 (10 probes). The predominant feature in both figures is the large horizontal scale of motion of the salt water wedge during each tidal cycle. If we use the 15 millimho/cm contour to separate the salt and fresh water masses, there are long time periods when salt water is not present at each instrument. This is certainly true at stations 2 and 3, although it is difficult to acertain at station 1 due to the large probe line deflections in February (Figure 13) and the loss of probes in March (Figure 14). However, the flood tide records at station 1 on both figures indicate a return of salt water. On flood tides probe line deflections were not large. No direct measurements of the salt water penetration length were made in this study; however, the stratification detected at station 3 on February 10 certainly indicates an intrusion distance exceeding 15 kilometers. The vertical spacing of the contours is indicative of the degree of mixing present and during most of the measured cycles the stratification is comparable on both the ebb and flood portions of the tide. This is only true for stations 1 and 2 in the March data; the results from station 3 indicate that the profiler was located close to the upstream limit of the penetration. The steepness of contours following the ebb tide show that the salt water is removed very rapidly and toward the end of ebb appreciable mixing of fresh water down into the salt layer has taken place. 47 The conductivity records in Figure 14 show that salt water did not move much past station 3 on each flood and as a result, the differences in stratification present at each profiler reflect the variation in mixing along the interface. The contour spacing generally increases from station 1 to 3 and suggests that both the degree of completion and rate cf mixing increase toward the wedge toe. If the mixing processes are predominantly two-way, this trend is reasonable since the salt water furthest upstream would suffer the greatest dilution and in turn become more susceptible to further mixing. Current Structure Simultaneous current metering was carried out at two stations on March 30, 1973. Profiles made 90 minutes apart (except for no. 6) are plotted in Figure 15 for stations 1 and 2. The times are indicated relative to the March 30 tide in Figure 15(a). The current speeds are plotted directly from the field data for station 1 (Figure 15(b)) since the flow directions were not measured. At station 2 the component speeds relative to the surface current are plotted . Throughout the period occupied at station 2 this surface current was ebbing and its direction varied only ±15° from the mean value. In general the current speeds are greater in both upstream and downstream directions at station 1 than those at station 2 due to the location of station 2 - on the inside of a sharp bend. An examination of the river bed topography near CL" < I < Q 0 STATION 1 k6 1 1 1 1 1 0 -150 -150 6 12 18 MARCH 30, 1973 (a) CURRENT SPEED ( CM SEC.) -100 -50 0 24 downstream upstream CURRENT SPEED (CM/SEC.) -100 -50 0 Figure 15. Current speed measurements made at tw stations in the Fraser estuary on March 30, 1973. 49 station 2 shows considerable sediment removal on the outside of the bend (Figure 4) and suggests the main flow stream to be located on that side rather than near the station 2 shore. One significant trend can be observed in the behaviour of the profiles on the return to ebb flows. After the surface currents have started to increase in ebb, the upstream transport of salt also continues to increase at station 1 (compare profiles 4 and 5). This has the effect of increasing the velocity shear as can be readily seen in profile 5. Such a trend is not so apparent in the station 2 data although the velocity shear also increases between profiles 5 and 6. These data also indicate a phase lag bewteen the tidal flows induced in each layer. Parameterization The conductivity charts provide a detailed description of the behaviour of the two water masses, but for further analysis and use in mathematical,modelling procedures, it is convenient to parameterize the data in terms of a length scale associated with the height of salt water. Since all the salt originally found in the water column came from the Strait of Georgia, we can decompose the measured water column into two layers: an upper layer of fresh water and a lower layer having an average salinity egual to that characteristic of the Strait of Georgia water. Thus the lower layer depth, h*, is defined by: h h» = VSg / S(z) dz (4) o 50 where Sg is the characteristic salinity of the Strait of Georgia water, S (z) is the measured salinity at depth z and h is the total depth of water at the probe line. In this way all of the salt in the original water column is redistributed into the salt layer of thickness h*. There are two main difficulties involved in this kind of parameterization. First, the probe lines do not extend to the deepest part of the cross section and so the h1 defined this way gives a salt layer depth relative to the bottom at each station and not referred to datum. At least 80 percent of the total depth was monitored at station 1 and the parameterization here is most accurate, having an estimated maximum error of +50 centimeters due to the portion of the water column not included in the integration. A greater error is to be expected at the other stations due to the probe line location (Figure 4). The second problem arises from the loss of deep probes which plagued the March data at station 1. The error introduced in this way is worst when the salt water level is low on the probe line, and the most reliable estimates are obtained near maximum intrusion. The calculated values for h' at station 1 have been superimposed onto the contour chart for a portion of the March 18 tide in Figure 16 and despite the loss of probes 1 and 2 the parameterization does reproduce the essential intrusion features once the salt wedge has reached to probe line. The reference salinity of 30.5 %c was used in the final analyses since it represents the Strait of Georgia water below 10 meters adjacent to the river mouth. The contours are lines of constant conductivity, not salinity; h'from equotion ( 4 ) Sg =30.5 %o 0 conductivity probes vertical scale: \.Q n- 2.0 meters 8 10 12 March 18 14 16 18 20 Hours Figure 16. Comparison of h» calculated from equation (4) and the measured conductivity structure on March 18, 1973 at Station 1. 52 however, the conductivity and salinity values are nearly equivalent for the temperatures found in the estuary. Figures 17 and 18 show the two-layer parameterization for data from both the February and March measuring periods. These figures span the same time periods as the conductivity data contoured in Figures 13 and 14 and allow a comparison between the calculated h* and the detailed salinity structure. In order to evaluate h1, a value must be found for Sg and for data presented here a figure of 28 %> was chosen. This corresponds to an average value of the maximum salinity recorded on the deep probes at station 1 over several cycles. The difference of 2.5 %o between this figure and the value of 30.5 %Q mentioned previously results in a difference in the interfacial heights not exceeding 25 centimeters. The conductivity charts indicate that the salt wedge is flushed out of the estuary on large ebb tides whereas the • interfacial line in Figures 17 and 18 never goes completely to zero at station 1. An examination of the digital data revealed that following the ebb tide, the entire water column was composed of salinities ranging from 1 to 5 Hi , which probably reflect the mixing of remnant salt water throughout the estuary. Station 1 is the shallowest section and there are several upstream pockets in which traces of salt water would remain after the main wedge has retreated downstream. Thus, for values of z less than 1 1/2 to 2 meters the interface shown in Figures 17 and 18 is an artifact of the integration prodedure and does not represent the salt wedge position. Bottom o o g' o ____ cd ' If) _a> o •x> OJ (m o o. CO o o p cvi O g ' o CO " .—. if) o QJ £ o —• O . CJ o o STN.2 Bottom STN.3 Bottom 10 _L_ 12 _i oooo HRS. FEBRUARY 1973 Figure 17. Two-layer parameterization of the conductivity data obtained in February 1973. Bottom lines represent the river bed at the platform sites. Free surface lines are plotted from measured data at Stations 1 and 3. on STN.I Bottom STN.2 Bottom STN.3 Figure 18. Two-layer parameterization of the conductivity data obtained in March 1973. Bottom lines represent the river bed at the platform sites. Free surface lines are plotted from measured data at Stations 1 and 3. 55 Although there is uncertainty in the actual interfacial position, dividing the water column into two layers this way provides a good graphic description of the phase relationship between the motion in each layer. Both figures show that maximum penetration lags high water, particularly at station 2, where the time difference falls in the range of 2 to 3 hours, consistent with the conclusion from the salinity profiles. The maximum penetration is not easily determined since the salt water often forms a plateau just following high water. One feature does stand out. There seems to come a point in the ebb tide when the salt water can no longer maintain the penetration, and after which the salt wedge is rapidly removed. For example, if we consider the data on February 10 (Figure 17), the salt wedge is completely flushed past station 1 (indicated in the contour charts in Figures 13) within three to four hours in spite of the considerable penetration recorded at station 3. This implies that velocities on the order of 60 to 80 cm/sec are required in the salt water. 56 Chapter 4. THEORETICAL CONSIDERATIONS An examination of the prototype data has revealed several important features of the salt wedge behaviour. It appears that large flood tides often produce intrusion lengths exceeding 10 to 15 kilometers above Steveston, and that salt water can be forced out of the estuary on extremely low tides, even for relatively small freshwater discharges. The two water masses also tend to retain their identifying properties throughout the tidal period in the estuary above station 1, despite the mixing created by the strong currents. The data also illustrate the distinct phase relationship between the salt and fresh water motions (Figures 13 and 14). A mathematical model of the tidal estuary flows which utilizes these data, can provide a quantitative relation between such important parameters as the salt wedge position and currents and the forcing functions of tide and river discharge. If the measured density structure at the station furthest seaward is decomposed into two distinct layers, as described in the previous chapter, the resulting information provides boundary data for a two-layer mathematical model. Since the upper layer is assumed to be fresh river water, the model does not allow mixing of salt or heat across the interface. A similar parameterization at the remaining stations provides data to verify the predictions. The field data also indicate that fresh water mixes down into the salt layer as it flows upstream reducing the salinity and hence the density contrast between the two water masses. 57 This effect would be greatest uear the wedge "toe'1 and probably produces a well mixed zone where the distinct stratification is lost. At any rate, detailed observations are not available for the leading edge of the salt wedge and in the subsequent theoretical work the wedge toe is considered to be defined by an arbitrarily chosen depth of salt water equal to one meter. These assumptions lead to a mathematical model of the estuary flows which may be visualized as shown in Figure 19. When salt water has been forced out of the estuary or in the river above the front, the tidal hydaulics can be computed by the usual barotropic equations which are well documented in the literature (Dronkers (1964,1970) ). The stratified or baroclinic calculation traces the motion of the salt front and the two-layer flows downstream of this position. The two calculations must be compatible at the boundary located by the wedge toe. Since the aim of the modelling is to examine the stratified tidal motions the governing differential equations retain the time derivatives. The final equations are difficult to solve analytically and recourse to numerical integration procedures has been made. To further simplify the equations lateral variations in the estuary geometry and density structure have been neglected. seaward limit sail front barotropic computation fresh water discharge Q Figure 19. Visualization of salinity intrusion into the Fraser estuary. freshwater density p river mouth Figure 20. Notation used in mathematical modelling. 59 The Equations of Motion In the case of a two-layer model without mixing, the fluid motions in each layer are governed by the laws of mass and momentum conservation. The equations of motion are referred to a right-hand cartesian co-ordinate system with the origin located at the river mouth (Figure 20 ). Neglecting lateral variations the longitudinal momentum equation for an infinitesimal volume in either layer can be written as (Cameron and Pritchard (1963) ): du + u dn + w d u = -a 3p (5) dt dx dz dx where u and w are the instantaneous velocities in the x and z directions, p is the pressure and 0! is the specific volume ( ) . In equation (5) the molecular viscous stresses have been neglected since they are several orders of magnitude smaller than the turbulent stresses. The instantaneous velocity u is considered to be composed of two components: (a) a time mean velocity, u obtained by averaging over periods of sufficient duration to remove turbulent fluctuations. Thus, u represents the slowly varying velocity field over time intervals longer than the averaging period. (b) a velocity deviation, u* , arising from the turbulent fluctuations having time scales shorter than used for averaging in (a). Similar relations apply to w. Further, if we neglect turbulent density fluctuations compared with velocity fluctuations, the mass conservation 60 equations are: du + d_v - 0 (6) dx dz d u + jTw. = 0 (7dx dz for both the instantaneous and mean velocities, assuming the water to be incompressible. Equations (5) , (6) and (7) combine to give the time mean longitudinal equation of motion: du + u du_ + ¥ du = - a dp - d <u*u*> - d <u* w*> (8) at ax dz ax ax di where < > indicates the time averaging process. The kinematic Reynolds stresses arising from the turbulent momentum transfers are represented in equation (8) by terms of the form -<u*Uj>. Pritchard (1956) has concluded, on the basis of his study of the James River estuary, that the only significant stress will be -<u*w*>, and that equation (8) can be further simplified to (dropping the overbar): du + u _du + w du = - et d_P " _d_<u*w*> (9) dt dx dz dx dz In order to derive the final model, equation (9) can be integrated over the depth of each layer, where the mean velocities and densities are considered to be uniform. That is, du/dz=0 in each layer. Onder these assumptions the pressure p will be hydrostatic, and the final equations are: (a) upper layer du + u du = - CL C dp dz - _1_ dt dx rj J, dx Tj where _dp_ = p g d h dx dx <u*w*> -<u*w*> \ (10) z=h z= h • 61 (b) lower layer J2U at h' 2 + u' au« = -a^rap'dz -_L_J<U*W*> ax h' J "air h • z=h • -<u* w*> z=0 (11) where ap' = pq dh + (yO1- yO ) q ah ' ax ax ax Since the turbulent stress terms cannot be explicitly evaluated in this form, a substitution will be required in terms of other flow variables. In this analysis the free surface stress will be neglected, and letting Ti = interfacial stress = -</3u*w*> at z=h' Tb = bottom stress = -<£>'u*w*> at z=0, eguations (10) and (11) become: au + u au = -g ah -Tx at ax ax prj au' + u* a u' = -g X a_h_ - ge_a_h_2 + Ti -Tb dt dx dx dx yO'h' (12) (13) where X =/o/yo'and € =(/o'-p)/p\ Equations (12) and (13) together with the continuity equations for each layer, form the governing system of non-linear differential equations: au + u au = -g ah - T± dt dx dx PV + _a_ (U V ) = 0 at dx a u' + u' a u' = -gX a h - g € dh' + ( Tj -dt dx dx dx p* h (14) ah' + _a_(u«h») =0 at ax These eguations are classified as hyperbolic in that they possess real characteristics in the x-t plane, and therefore impose definite requirements for the presentation of boundary and 62 initial data. No general analytic solutions have been found for eguations (11) even for fluid flows where the turbulent stresses are negligible. If the non-linear terms are excluded the method developed by Rattray (1964) may be used; however, in this application there is no a priori reason to suspect the convective accelerations are negligible. Approximate solutions have been obtained in this study using explicit finite difference techniques, although other numerical procedures such as the method of characteristics (Abbott (1966) ), or implicit finite differencing (Grubert and Abbott (1972) ) can be used. Equations (14) form the basis for a mathematical model of stratified flows, but the appropriate boundary data and initial conditions are also required before a solution can be obtained. In addition the boundary stresses must be formulated in terms of the four dependent flow variables and these considerations are discussed in the subsequent sections. Thg Characteristic Structure and Boundary Conditions If the momentum equations are written in terms of V) and a It of h and 1 u -n u« + a ax h« u2/2 + g ( 77 +h») u i) u,2/2 + g(h» + Xi7 ) u'h« - Ti/y077 0 ( Ti- Tb) P' h' (15) 6 3 Or d_V + dF (V) = f at a x (16) where v, F (V) and f are vectors as defined in equation (15). Since V = V(x,t) the equations of variation are: (17) dV = _r3_v dt + aV dx at dx Written together in matrix form, equations (16) and (17) become: or ulgOOOgO V 0 u 1 0 0 0 0 0 0 gX 0 u' 1 g 0 0 0 0 0 h' 0 u' 0 dx dt 0 0 0 0 0 0 0 0 dx dt 0 0 0 0 0 0 0 0 dx dt 0 0 0 0 0 0 0 0 dx dt [A] [B] = [C] 2nx2n n n au/ dx du/ at dV/ dx dV/ at au1/ dx auv at ah'/ dx ah«/ at fl 0 f. 3 0 du AV du' dh' (18) where 2n is the number of dependent variables and n is the number of discrete flow layers. Abbott (1961,1966) has shown that the characteristic values for this system of equations are given by the determinant [A] = 0. After elementary operations on the matrix [A], g 0 g (u-c) 0 0 gX (u'-c) g (19) 0 h« (u«-c) c =x=dx/dt, the inverse slopes of the characteristic curves in the x-t plane. I A | = (u-c) V 0 0 64 Finally: [ (u-c)2 - g 77 ][ (u •-<;) 2 - gh' ] - g2 X 77 h' = 0 (20) Equation (20) is the characteristic quartic (for a two-layer system) whose four roots divide into pairs associated with each layer. Physically each pair of roots represents the speed of long gravity waves on their respective surfaces. These roots can be evaluated numerically, or approximated to order e by the following relations (Schijf and Schbnfeld (1953)): c1 = u + ± / gh (21) h J cf = uh* + u'V ± / ge 77 h' - (u-uM 2 rj h' (22) h • ' h 2 where c* represents the positive and negative surface gravity waves and cf , the positive and negative internal waves. Now, in two-layer flows described by four dependent variables, four boundary conditions must be known for all time. Furthermore, the boundary conditions must be supplied to the domain of computation in accordance with the characteristic structure. To guote Grubert and Abbott (1972): "The boundary data is of m-point type where m is the number of characteristics entering the region of computation from the boundary in the direction of computation." Three situations are of interest for the problem under consideration and are summarized in Figure 21. In each graph the two-layer region of computation in the x-t plane is shown in part of which the equations have been solved. We may imagine that the solution is being advanced in time to the two points P, and P2 65 Subcritical flows at both boundaries .'. two boundary conditions required at each end. x=0 x = toe Solution determined Supercritical outflow at both boundaries I boundary condition at seaward end,and 3 boundary conditions at upstream end. x=0 x = toe x= 0 Solution determined x = toe Supercritical inflow at both boundaries I boundary condition at upstream end , and 3 boundary conditions at seaward end. Figure 21. Possible flow states for the two-layer stratified computation. The required number of boundary conditions is indicated. 66 using the method of characteristics (Smith (1965) or Abbott (1966) ). For each point the four characteristic curves have been sketched corresponding to subcritical flow, supercritical outflow and supercritical inflow. Since field observations show that tidal conditions in the estuary do not result in a surface bore, (that is, c+>0, and c~<0 for all flows) these flow conditions refer to the propagation of internal waves. The specification of the boundary values is shown in Figure 21 and applies regardless of the solution technique. A Steady Flow Model - Initial Conditions Any solution of equations (14) requires initial or starting values for the dependent variables. In the numerical procedures used . by Vreugdenhil (1970) and Boulot et al. (1967) the unsteady equations of motion were transformed to the equivalent steady flow relations and then numerically integrated to provide the initial conditions. An alternate approach, adopted here, is to make use of an analytic solution for the "arrested salt wedge" to provide the required values. This method provides insight into the mechanisms controlling, or allowing the existence of a stationary salt wedge in suitable estuaries. Steady-state models for salt water intrusion have been studied in the past and the mechanisms governing the wedge behaviour described, principally by G. H. Keulegan (1966). In his analysis Keulegan deals mainly with the penetration length of 67 arrested saline wedges and evaluates the necessary forces from experimental data obtained in laboratory flumes. He arrives at the following equation for the intrusion length in an estuary: L0 = 6.0 h R! (2FA )* (23) where LD= intrusion length, h = total water depth, RA = densimetric Reynold's number = J g e hyC , v - kinematic viscosity of water, FA = densimetric Froude number = U (L0)// g e h U (L0) = fresh water velocity at x=L . Equation (23) was derived from dimensional considerations and experimental data. In this way no assumptions were required regarding the interfacial stresses which provide the arresting forces on the wedge. Keulegan found the shape of arrested wedges to be similar over a wide range of river parameters and to be a function of x/Lo, although this relation is not given. A second analysis of an arrested salt wedge was described by Farmer and Morgan (1953). They present a closed form solution for the shape of the.arrested wedge by assuming the interfacial stress varies as K p 0 (L)2 and will remain constant along the wedge. To evaluate the solution K must be determined from the penetration length, as in flume studies, or from the wedge geometry in prototype situations. In effect, then, the solution is '•calibrated" to the data using the coefficient K. It is also possible to derive an analytic solution for 68 the wedge shape in the following way. If we assume the salt wedge to be stationary and that no mixing of salt and fresh water takes place, the governing equations (14) reduce to: udu + gdh + Ti = 0 (24) dx dx pi) dJUT/ ) = 0 (25dx g X dh + g e dh' - Ti = 0 (26) dx dx P'h* with the notation as defined in Figure 22. Either a substitution will be required for ri or it can be eliminated from the equations. Following the latter method equations (24) and (26) reduce to: h - h' d (uZ) + q d (h?) +^e_d_(h»2) = Q 2 dx 2dx 2 dx where €={p'-p)/p' . Letting V - h - h* and making use of the fact that steady flow exists in the upper layer, that is a constant discharge Q = u 17 (from equation (25) ) we have: Q2_d(!/7?) + g_d_(hZ) + _g_e_d_(h •2) = 0 dx 2dx 2 dx where each term is an exact differential with respect to x and thus integrable. The final result is: Q2 + gh£ + q 6 h '2 = C (27) 77 2 2 where C = constant of integration. Solving for h»: hi3 - hh'2 + (h2 -C)h' + (Ch - h3 - 2Q2) = 0 (28) e e g ^ The constant C is determined from the boundary conditions, either at x = 0 or at x = L where h'= 0. Now, in general the distance L is unknown and at x = 0 the depth of salt water is 69 Fresh water density /O z = h x = 0 x= L, Figuce 22. Notation used in steady flow solution. 0.75 to O xrixO.5 0 0.2 5 0 o Barotropic solution K = 0.006 Discharge Q (m/sec.) Figure 23. Stags-discharge relationship Fraser estuary. for 70 unkown. However, if a control section is postulated for x = 0, the equation for the characteristic wave speeds provides the required relation. This appears to be a reasonable hypothesis on the basis of Stommel and Farmer's (1952) investigation showing that sudden increases in channel width can produce control sections in steady two-layer open channel flows. The concept of a control section comes from flow conditions which prevent progressive waves from propagating past a certain point; waves moving against the current become stationary just at the control. In the case of two-layered flows an "internal" control section prevents the upstream propagation of long internal waves and in effect means that no communication from flow conditions in the wider receiving waters to those in the narrower channel can occur along the interface. Only the positive surface gravity wave can pass the control in the +x direction. Thus equation (20), with u' = 0 and c| = 0 becomes: (uz - gh) (-gh) - g^XT; h' = 0 where u = Q/y , .-.gey - °-2/ V 2 h« = h - (Q2/g e )'/3 (29) The constant of integration C can be determined from conditions at x = 0, where if h'(x=0) = a (Figure 22): C = €a2 + h2 + 2Q2/g(h-a) (30) Equations (28) (29) and (30) enable us to specify the salt wedge shape if the surface elevation is known; indeed the only x - variation in equation (28) enters through the surface term h. This derivation is substantially the same as described 71 by Hodgins and Quick (1972) although it has been simplified by the choice of the bottom as level datura and some new computations will be presented. In this previous study equation (28) was given in non-dimensional form by dividing with a dimension equivalent to h (L) : R3 - R2 + (1 - Ci )R + (Ci - 1 - 2FA2) =0 (31) where R = h/h (L) C, = C/h (L)z = densimetric Froude number = u (L)/y g e h (L) It is interesting that FA appears directly as a result of the convective acceleration term, and does seem to be a useful and fundamental parameter of stratified flows. In most unsteady tidal computations the convective acceleration term, udu/dx, can be neglected as it is usually one or two orders of magnitude less than du/dt and therefore we might expect a similar situation for steady flows. This hypothesis can be tested by solving equations (24) and (26) without the udu/dx term, in which case the solution for h' has the form: h2 + € h'2 = C (32) an elliptical equation. As before C can be evaluated at x = 0 with the use of equation (29): C = e a2 + hz (33) solving equation (3 2) for h': h' = +/ (C - h2)/€ (34) which has a real root only if C > h2. Thus the value of x at 72 which h2 = C is the "penetration length" of the salt wedge in this very simple analysis. It must be noted that h and Q are not independent. Together they represent the unique stage - discharge relationship for a given estuary. To specify h and Q and in turn evaluate h', I have used a barotropic model derived and solved in the following way. If the momentum eguation (9) is integrated over the total depth and coupled with the continuity equation a barotropic hydraulic model is obtained: au + u _£a = -g _3JL - rb (35) at ax ax ph dh + _a_(uh) = 0 (36at ax Eguation (36) is a continuity calculation on a per unit width basis and provides only a crude approximation to rivers having irregular channel forms. To improve the simulation the continuity equation is often applied to the whole cross-sectional area, A, in the form: b dh, + _a_(uA) = 0, (37) at ax where b = storage width. The bottom friction term ^b was replaced with the relation Kp u|u| where K is a variable dimensionless friction coefficient, and then the equations were calibrated to the tide and discharge conditions for six days during the March measuring period. The equations were solved using an explicit central difference technigue similar to the scheme described by Cronkers (1970). The unit width continuity equation was used here for ease and 73 speed in computing; details of a more sophistocated computational procedure based on equations (35) and (37) appear in following sections. Once the friction coefficients were obtained the model was run for a zero amplitude tide and varying discharges to produce the mean surface slope in the lower estuary. The relation between Q and dh/dx obtained in this way is shown in Fiyure 23. Compared with the measured tides at stations 3 and 4, the calibrated unit width equations were successful in predicting the surface elevations within ± 5 percent at low water (the worst fit and corresponds to a difference of approximately 70 cm) and within ± 20 minutes for high and low water times (Figure 24). The two most important parameters are dh/dx and the density contrast e . In Figure 25 (a) four typical solutions of equation (28) are shown for a variation in Q and hence dh/dx, with € =0.025, while in Figure 25(b) the parameter € is varied for constant dh/dx. Solutions of equation (34) are also plotted in Figure 25 (b). The penetration lengths are increased with an increase in € and a decrease in dh/dx. The solutions excluding u du/dx have shorter penetration lengths than those incorporating the convective acceleration since this term resists the freshwater velocity changes in the upper layer, having the same sign as the interfacial stress. Thus the same discharge in each solution can be obtained with a lower stress when udu/dx is retained. The interfacial stress arises from the growth and breaking of internal waves associated with the region of large density gradient, a process enhanced by increasing velocity shears. Ultimately the convective acceleration reduces the 2r Station I x = OKm Station 4 x = 24.4Km / / + 12 _L_ V G.D G.D -\T— G. D 18 _l_ 0 12 1.8 0 Hours March 17 • o + Measured values from tide gauges March 18 Barotropic prediction Figure 24. Measured and predicted surface elevations at three stations on the Main Arm of Fraser estuary. 75 20 30 x (Kilometers ) (a) dh xio6 40 50 20 30 x ( Kilometers ) 40 50 (b) Figure 25. Typical solutions for an arrested salt wedge. 76 velocity gradient and hence reduces the stress compared with a hypothetical system without such an acceleration term. Once the interface is established the stress, Tx, can be calculated from eguation (26). Figure 26 shows the variation of ri along the wedge increasing from near zero to a maximum toward the upstream end. The stress then rapidly decreases to zero at the toe. It must be pointed out that these stresses are calculated from a uniform surface slope and may not accurately depict the situation in real estuaries. Since the velocity shear would be increasing toward the mouth an increase in Ti is expected. This would result in an increasing surface slope to provide the necessary pressure force driving the upper layer flow. The decreasing upper layer thickness and increasing velocity produce the right conditions to form a control section near the river mouth. The stresses calculated here serve to fix an order of magnitude to these forces, at least in the vicinity of the wedge toe. Farmer and Morgan (1953) assumed the stress varied as K^ou(L)2 and using experimental or prototype data, evaluated K for particular examples. They briefly discuss fitting their wedge model to the Mississippi River and state that a K of 0.001 was appropriate. Using their data I find this corresponds to an average stress of 1.8 dynes/cm2, which compares favourably to the larger values shown in Figure 26. The South West Pass of the Mississippi River has a densimetric Froude number of approximately 0;25, derived both from the data of Keulegan and 77 x/L Figure 26. Typical interfacial stress curves for arrested salt wedges assuming a uniform surface slope. 78 Farmer and Morgan. Vreugdenhil (1970) used a K of 0.005 in his unsteady analysis of the Rotterdam Waterway, which compares more closely with the value Farmer and Morgan report for their experimental studies (K = 0.006). The larger value might be expected in the unsteady case where turbulent salt water flows would enhance mixing between the water masses. The stress formulations are not identical since Vreugdenhil allows the stress to vary as a function of the local velocity shear squared whereas Farmer and Morgan assume it is constant along the wedge and a function of the freshwater velocity squared at x=L only. In summary it is possible to solve the steady flow equations by eliminating the interfacial stress terms, but in the absence of information about variations in the surface slope, final conclusions regarding the form or variation of the stress cannot be derived. A uniform slope gives maximum stresses of about 2 dynes/cm?, in agreement with other studies, but provides a longitudinal stress variation conflicting with the concept of increased mixing due to the formation of a control section near the river mouth. The importance to this study lies in equation (28) which provides the required initial values for solving the unsteady model. It now remains to formulate the boundary stresses Ti and Tb in terms of the dependent flow variables. Boundary Stresses The Reynold's stresses resulting from the turbulence in each layer appear in equations (14) as the bottom and interfacial 79 shear stresses and represent the closure problem associated with all turbulent flows. In this case there are four equations in six unknown quantities, since there is no explicit way of representing -<pu*w*> at the interface or bottom. Traditionally, an assumed dependence of -< pu*v*> on the other dependent flow variables is substituted in the equations of motion. The bottom stress, ^b, is often assumed to be like a drag force which varies only with the mean layer velocity squared. Both Vreugdenhil (1970) and Boulot et al. (1967) replaced Tb by relations of the form K yO'u ' | u ' |, K being a dimensionless coefficient of proportionality. This formulation assumes that changes in the rate of turbulent momentum diffusion are adequately represented by the u'2 variation and neglects any length scales associated with the boundary generated turbulence in the lower layer. I have also used this relation for Tb, both in the barotropic and stratified equations. Turbulent mixing in the presence of a density gradient is not a well understood subject due to the complexity of mechanisms which may be involved and the difficulty in generalizing between specific examples. The Reynold's stresses arise from the vertical turbulent diffusion of horizontal momentum and depend upon the intensity of turbulence on either side of the density interface. Considering the two layer system, if velocities in the lower layer are sufficiently large, turbulence generated at the bottom may have large enough scales to produce mixing at the 80 interface. When the scales of the wall turbulence are not on the order of the layer thickness, viscous dissipation considerably reduces this effect. Turbulence may also be generated at the interface itself due to the velocity shear. In this case the turbulence is associated with the breakdown of internal waves or the growth of Kelvin-Helmholtz instabilities (Turner (1973)). Biles (1961) proposed that away from the influence of boundaries a sufficient condition for stability to small disturbances in a inviscid, continuously stratified fluid is that the gradient Richardson number, Ri = N(z)/( dn/ dz)2 > 1/4 everywhere in the flow. N (z) is the Brunt-Vaisala frequency, given by -q(dp/ dz)/p , and approximates the frequency of oscillation of a water parcel given a small initial displacement in a stable water column. Turner (1973, $4.1) has pointed out that this does not mean that instability always results when Ri falls below 1/4 and has referred to some examples. Alternatively, it appears that turbulent mixing can persist even for Ri » 1/4, as shown in the results discussed by Taylor (1931). In analyzing continuous velocity and density profiles from a tidal channel, he reported eddy viscosity coefficients, defined by V€= -<uV*>/(du/dz)f on the order of 1.0 to 3.0 cm2/second for Richardson numbers of ten or more. Compared with a kinematic molecular viscosity of order 0.02 cm2/second, it seems that turbulent momentum diffusion can persist for Richardson numbers exceeding 1.0. It is interesting to compare some of the Fraser River observations with those of Taylor and the critical Ri of 0.25. I have taken the velocity 81 profiles from March 29 at station 2 and computed the velocity shear, du/ dz, using a Az of one meter. Since each profile required about 10 minutes to complete, the gradient Richardson numbers are "time averaged". Furthermore the density determinations were not made until the completion of velocity profiling. The Brunt-Vaisal'a frequency, N (z) , was also computed using a Az of one meter. The results are listed in Table 3 and plotted in Figure 27. The pycnocline was quite diffuse on both observations but the velocity shears were sufficient to produce Richardson numbers on the order of 0.5 to 1.0, and in one case, less than 0.25. Values of Ri on this order and the comparison with those of Taylor indicate that turbulent mixing generated internally, is to be expected in the Fraser estuary. Although velocity and conductivity data were collected at two stations on March 30, I have not attempted to calculate eddy viscosity coefficients from the momentum equation since the probe lines had suffered considerable damage and the density was specified at too few points. Furthermore, in a river with as many bends as the Fraser, secondary currents are strong enough to make even the straight salt balance based on two stations of doubtful validity. As far as I know direct measurements of the Reynold's stresses, and determination of their relationship to velocity and density profiles have not been made in stratified estuary flows. Such investigations have been confined to laboratory experiments where the usual practice has been to calculate -p <u*w*> by measuring the remaining terms in the equation of motion. Ellison and Turner (1959) found the rate of entrainment by surface jets 82 Table 3 Comparison of Richardson numbers obtained from Fraser River data and Taylor (1931). N (z) =-gA/o / p {\z is the Brunt- Vaisala frequency, Ri=N(z)/( Au/ Az)2 the Richardson number and V€ is the kinematic eddy viscosity. Depth Au/ Az N(z) Ri (meters) (sec-1) (sec-1) (cm2/sec) March 29, 1973 (0830 hours) station 2 1.5 0.10 0.0078 0.78 2.5 0.30 0.0176 0.20 3.5 0.26 0.0578 0.85 4.5 0.23 0.0392 0.74 5.5 0.08 0.0166 2.77 March 29, 1973 (1552 hours) station 2 0.5 0.03 0.0039 4. 33 1.5 0*10 0.0127 0.75 2.5 0. 25 0.0235 0.38 3.5 0.18 0.0461 1.42 4.5 0.23 0.0235 0.44 5.5 0.09 0.0539 6.65 Shultz• s Grund Taylor (1931) 2.5 -0.010. 7.4x10-6 7.14 3.1 5.0 -0.017 11.0 3.85 3.1 7.5 -0.022 27.9 5.88 2.7 10.0 -0.024 58.9 10.2 2.2 12.5 -0.019 103.0 28.6 1.9 15.0 -0.008 80.8 125.0 3.8 17.5 0 45.6 83 crt =(/D_ |.0) x I03 U (cm/sec) .5 10 15 20 25 30 0 20 40 60 80 100 120 u> 2 CD CD sz Q_ , CD 51 Q 0830 March 29 Station 2 Az= I m. *>l.2 10 20 30 40 5 0 60 0 0.2 0.4 0.6 0.8 1.0 12 2 N(z)=-x I03 (Rad/sec) Ri(z)=N(z)/(AUy 10 15 20 U (cm/sec) 25 30 0 20 40 60 80 100 120 10 20 30 40 50 60 0 0.2 0.4 0.6 0.8 1.0 1.2 N(Z) Ri(z) Figure 27. Profiles of °"t and velocity on March 29, 1973 at station 2. Distributions of the Brunt-Vaisala frequency and gradient Richardson number are shown. 84 and inclined plumes could be related to an overall Richardson number, RiQ = g e h/U2, where h is the layer thickness and U is the average layer velocity. The entrainment coefficient, E, equal to 1/U [d(Uh)/dx], was determined in two series of experiments and found to be non-linearly related to Ri0, falling off sharply as RiQ approached 0.2 to 0.3. Entrainment was effectively zero for Ri0 > 0.8 to 0.9. A relationship between the Reynold's stress at the interface and the other • flow variables does not emerge from this work; however, a dimensionless number such as Ri0_1 or Ri0—'^2 would appear to characterize the mixing properties of the flow. In a series of flume experiments with a turbulent layer of salt water flowing beneath still fresh water Lofquist (1960) showed the shear stress at the interface depended upon both an overall Richardson number and a Reynold's number, Re, associated with the salt layer. The inverse square root of Ri0 is the interfacial Froude number, Fi=U/ J g • h, g'= eg (Stommel and Farmer (1951)) and ft was observed to increase with increasing Fi and Re. Lofguist's study was primarily aimed at relating -<u*w*> to the velocity and density profiles in the flow and he did not attempt to derive a relation between the stress and the integrated layer properties. The laboratory experiments have teen confined to steady flows over (or under) an ambient fluid of almost negligible velocity, and turbulence would be well developed only in the plume. Since stratified flows in the Fraser estuary are unsteady and both layers can have significant velocities, extrapolating the experimental results is neither possible or reasonable. 85 However an interfacial Froude number does seem to characterize the flows and parameterizing the stress, ri, in terms of it has an observational basis. The interfacial Froude number can be interpreted in terms of the internal waves on the interface. If the ratio of h */h « 1 (as in fjords, where u' ~ 0), J 9'V (Figure 20) is the approximate phase speed of infinitesimal internal gravity waves on the density discontinuity in the absence of currents. Then, Fi=0/yg, where U = u-u' is the relative upper layer velocity, compares the phase speed of internal waves and the layer velocity, for if Fi=1, the internal waves cannot travel against the current and vigorous mixing ensues. When h'/h<1, the phase speed estimate can be improved by using J q%rq h*/h, and writing Fi=U/|// g»i7 h*/h. The Froude number defined this way is not meant to accurately depict the conditions under which the internal waves will be arrested; as shown in equation (22) the phase speed is a complicated function of the layer depths and velocities. Instead Fi possesses the correct functional form to simulate the approach of critical conditions. Both Vreugdenhil (1970) and Boulot et al. (1967) have assumed the interfacial stress is proportional to Kp U|U|, an expression which has the dimensions of stress and allows for a dimensionless friction coefficient. Such a formulation includes only the relative layer velocity, u-u* and makes no reference to changes in layer thickness which also influence the speed of internal waves. The Fraser observations show that the fresh water depth can become very small near the estuary mouth and that 86 the salt layer thickness decreases near the toe. These variations, which reduce the interfacial wave speeds and promote mixing, can be included in the stress formulation using Fi. Since one of our study objectives was to examine the influence of various stress functions on the form of the predicted salt wedge shape and the residence time of salt water in the estuary, four substitutions for Ti have been made: (i) Ti = Ki p Fi (38) (ii) Ti = Ki p U J U | (39(iii) Ti = Ri pUIUIFi (40) (iv) Ti = Ki ^UIDIFi2 (41where each formula preserves Ki as a dimensionless coefficent except (i). In a no-mixing two-layer model the denominator of Fi becomes what may be called a weak function of the layer depths having a maximum when - h' = h/2. In the numerical solution V and h1 never have the limiting values of zero. Thus Fi takes on the role of an amplification factor making the stress more sensitive to the velocity shear and increasing when either layer thickness decreases rapidly. Relation (i) has been included since Ti is often assumed proportional to u-u' in modelling stratified flows. 87 Chapter 5. SOLVING THE TWO-LAYER MODEL By recasting the stratified flow eguations (14) in explicit finite difference form, and combining them with a barotropic model of the entire estuary, solutions were obtained for the tidal motions of the salt wedge during the observation periods. The accuracy with which a numerical model simulates the prototype behaviour is derived principally from the accuracy in specifying the boundary conditions. Observations of stratified flows show that for densities similar to those of sea water, the free surface barely feels the presence of the baroclinic flows. For example, in a series of laboratory experiments on gravity currents in which a barrier dividing a long flume into, two channels of fresh and salt water was suddenly removed allowing the heavier fluid to flow under the fresh water, the free surface was immeasurably disturbed by the current's progress despite occupying about one half the water depth. In the estuary it is reasonable to expect the free surface to be nearly independent of the salt water flows, driven only by the tide at the Strait of Georgia and the upstream discharge. Under this assumption the barotropic model provides both the free surface and the net discharge at any section and the baroclinic flows can be computed considering only the lower layer eguations. This also means that the calculations are perfectly matched across the salt front since the free surface and fresh water discharge on either side come from the barotropic model. This chapter presents the finite difference equations and a brief description of the barotropic model, along with a discussion of the computational procedure and 88 application of the boundary conditions. The computed salt wedge motions are then compared with observations and some essential features of the model are examined. The Finite Difference Eguations Approximate methods for solving hyperbolic equations fall into three broad categories: the method of characteristics which utilizes natural properties of the differential equations and explicit or implicit finite difference schemes. There are no dominant reasons for choosing between the finite differencing techniques, both have strengths and weaknesses. Explicit central differencing has found repeated application in solving many kinds of wave equations, while Vreugdenhil (1970) and Boulot et al. (1967) have made use of "Lax-Wendroff" techniques (Richtmyer and Morton (1967) ), utilizing forward differences in time and central differences in the spatial variable to solve the two-layer flow problem. Explicit schemes are attractive from a programming point of view but suffer from restrictive stability conditions. Grubert and Abbott (1972) have successfully applied an implicit scheme to the stratified flow equations which is not confined to small time increments by stability considerations. In general, however, implicit schemes require a more involved algorithmic structure. The one-dimensional hydraulic model using equations (35) and (37) in central difference form was developed for the entire Fraser estuary including all major branches and Pitt Lake. In 89 the discussion to follow i will reference a time line in the x-t plane and j will refer to a space or x position. Consider first the continuity equation (37) in central difference form: bj (h i + i -h i-l 1 > 2 St j-t Ajll -U j-l Ai-2 Sx (42) Since Ajj|,Aji| are unknown they are replaced by the corresponding values, obtained by arithmetic averaging, from the preceding time row; that is: )/2 = A);1, A 1  Aj+I (A]^+A;-A ' — (Ij-jtA1-1 )/2 = Vr\ the overbar here denoting the average. Thus: hj|' = h1-:1 - 1 St J bj Sx (43) The computation molecule is illustrated in Figure 28 (a) and shows exactly what information is known and utilized for advancing equation (43) •. Proceeding in a similar manner for the momentum equation (35) (Figure 28 (b) ) , where the bottom stress rb= Kb/0 u1}1 lu1]1 | : u-'j1 =< i-l 1 + |i.(uj^ -ujj>) • 2 StKb lu'J1 | /hj 2 Sx Equations (43) and (44) comprise the barotropic model. (44) 90 initial conditions J"2 J"' J j+l j+2 (a) advancing the continuity equations it initial conditions - x J*» J + 2 (b) advancing The momentum equations h ,h'known ©u,u'known Qh.h' to be determined rju,u' to be determined Figure 28. Finite molecules. difference computation 91 For simplicity the stratified computations are done on a unit width basis and the lower layer continuity equation has been written in a difference form analagous to equation (43): h' ' + ' = h» !T J j u'.' . h -u« .' , h» I-1 J+I J+I j-l j-i (45) with the computation molecule in Figure 28 (a). Using equation (40) for the interfacial stress and the quadratic friction law for the bottom stress, the difference momentum equation is: u'1'}1 =iu»iJl - St 8x I I IT i-l X(hj|, -!»._!,)+ eCh-jj, -h'.j., ) 28t(Ki U'T' | u |Fi - Kb u» '-' | u«! j 1 |/h» j )l/ S t uij-1 -u'i-' + 1 where U '. 1 = (u-u*)^1 = velocity shear, and Fi (46) •i = oV // ge V j h» j /h-and is computed on the molecule in Figure 28(b). If the two continuity equations from (14) are summed: (ui? ) + (u»h') d (V +h«) = _a_h = - a at at dx Since h is identical in stratified and barotropic computations, dh = - - a (.b.h) - - (UT/ ) + (u'h«) at ax ^ " / ax where ub = barotropic velocity. Thus solving for u't1 : i4lri (47) the constant of integration being zero. The five equations (43) 92 to (U7) form the basis of the stratified flow algorithm and as can be seen from the computation molecules the minimum number of grid lines required for a computation is six, which serves to define the minimum penetration length that can be admitted. This is a limitation arising from the central differences and is unfortunate since reasonable values for Sx are on the order of one kilometer and the separation of stations 1 and 2 is only 2.7 kilometers. Thus any error in establishing the initial conditions is immediately felt at station 2. Explicit difference schemes for hyperbolic equations are subject to restrictions in the choice of the grid spacings Sx and St , which are directly related to the characteristic structure. For example, in Figure 29 the shaded curvilinear triangle PQR defines the continuum domain of dependence for the solution at point P, and a unique solution is determined inside this domain (Crandall (1956) ). The curves PQ and PR are the dominant characteristics, the c+ and c- characteristics for the stratified flow equations. The fundamental paper on the numerical treatment of hyperbolic equations dates back to 1928, by Courant, Friedrichs and Lawy, in which they state what has come to be known as the Courant-Friedrichs- Lewy (CFL) stability condition: in a numerical grid which does not follow the characteristic directions, the finite difference domain of dependence must at least include the continuum domain of dependence. Stated in other terms, the slope of the finite difference characteristics must be less than or equal to the slope of the characteristic curves. In Figure 29 this requires Figure 29. Illustration of relationship between continuum and finite difference domains of dependence for a stable explicit difference scheme. 94 that the finite difference solution point S, advanced on information from points Q and R, must lie inside the continuum domain (below P). This stability condition has a simple physical explanation: we must integrate the numerical equations faster than the maximum celerity of gravity waves in the flow. Otherwise information used to advance the solution would be out of date before the spatial integration is completed within each pass. Quantitatively, these restrictions are ensured if St < Sx / max (c+) (48) where max(c+) is the largest value of the surface gravity wave which will occur. In this study St has bean chosed as close to the value of Sx /max (c+) as practicable. Rigorous mathematical theories for convergence and stability are' well developed only for linear systems of equations (Roach (1972)). This means that we are dependent on experiment and intuition to assess the accuracy of numerical solutions to systems such as equations (43) to (46), and serves to underscore the need for field data. A convergent finite difference scheme is defined as one whose solution approaches the exact solution of the differential equation as Sx and St approach zero. Since the difference equations obtain the solution through a series of arithmetic operations of finite accuracy, round-off errors are introduced into the dependent variables at each point. Stability is usually defined in terms of the growth or decay of these errors (Smith (1965)). In a dissipative and stable difference 95 scheme, numerical errors either decay exponentially or are bounded to sufficiently small values. Unstable difference schemes are always non-convergent. Due to the difficulties in proving convergence and stability for linear systems with constant coefficients, I have not attempted to derive convergence criteria for the non-linear equations involved in this study. Instead, the accuracy of the solutions is judged in terms of the simulation of the field data within the assumptions made in deriving the equations. The Barotropic Model The one-dimensional or barotropic model evaluates equations (43) and (44), starting at the various river mouths and integrating upstream to the end of tidal influence. The estuary has been schematized into seven major branches partitioned into one kilometer segments as shown in Figure 30. The heavy outlines indicate these branches while the lighter lines show the remaining channels not included in the scheme. The most notable omission is Sea Reach between branches (1) and (3) which does have an appreciable storage. The cross-sectional areas, conveyance widths and channel depths were obtained at every grid point except those in Pitt Lake, from the sounding data provided by the Department of Public Works (Government of Canada). In practice the conveyance width was computed over the portion of river having depths greater than two meters at low water and the channel depth was calculated at 85 percent of the greatest sounding. Detailed bathymetric data are not available for Pitt 96 Figure 30. used in the Schematization of Lower Fraser barotropic hydraulic model. R iver 97 Lake, so a shore to shore conveyance width was taken from maps, but depths and areas were estimated. Since §x equals one kilometer St < 83 seconds by equation (48) and for convenience in the programme bookkeeping, a St of 60 seconds has been used. The boundary conditions applied to this model were the measured surface elevations at station 1 and the Agassiz discharge at the final velocity section of the Main Arm. In the absence of detailed discharge data at various points along the estuary (from which the friction coefficient Kb can be determined), the model was calibrated to measured surface elevations at stations 3, 4, 5 and 6. A similar model based on spatial increments of 5000 feet was extensively calibrated against the data of Baines (1952,1953) by Joy (1974) in his study of dispersion in the Fraser River. To provide additional comparisons the predictions from Joy's calibration were also incorporated. The predicted and measured surface positions for March 30, 1973 are plotted in Figures 31 and 32. Within the lower Main Arm (below station 4) the phase agreement is satisfactory, however there is a significant error just after the turn to ebb following higher-high water at station 3. I attribute this to the crude simulation of the storage capacity along the south bank of the river below station 3 and the influence of the junction to Canoe passage. Both the phase and elevations are accurately predicted at station 4 (within 10 minutes and 5 centimeters of elevation). The predictions appear to lead the measure response at station 5 on the order to 30 to 40 minutes, suggesting that the depths 98 1.00 r -I.OOL Station 3 Barotropic prediction Sx - I Kilometer • Sx = 5000feet o observed G.D. G.D. 8 10 12 14 16 March 30,1973 18 20 22 24 Hour, Figure 31, Comparison of predicted surface elevations with obervations at stations 3 and 4. 99 i.oor 0.75 Station 5 Barotropic prediction ~ 8x = I Kilometer • Bx- 5000 feet o observed G.D. G.D. 8 10 12 14 16 March 30, 1973 18 20 22 24 Hour, Piyuce 32. Comparison of predicted surface elevations with obervations at stations 5 and 6. 100 downstream are too large and that the complicated storage areas (the river is very wide with extensive banks covered only at high water) and the influence of the Pitt River junctions are poorly modelled. The phase relationship - is satisfactory at station 6,; however, the high and low waters are over-estimated by approximately 10 centimeters. The velocities measured at station 1 and 2 were spatially averaged and compared to the predicted velocities at sections 3 and 5 (Figure 33). Although the measurements are very limited, the phase agreement seems satisfactory but the magnitude of large ebb currents is underestimated, especially at station 1. This may be due to neglecting the storage in sea Reach which would ebb into the area of Station 1 or an overestimate of the conveyance areas. I have calibrated this model by altering only the friction coefficients although Joy (1974) found that changes in the channel geometry were also necessary to achieve the best fit possible within the limitations of the schematization. As a tool for investigating the stratified flows, the barotropic model calibrated as in Figures 31 to 33 was satisfactory. Computational Control and Boundary Conditions Once the barotropic model was calibrated, the stratified flows were computed by interfacing the lower layer difference eguations with the basic model. All finite difference eguations were coded in Fortran IV and solved using the IBM 370 Model 168 computer available through the University of British Columbia 101 Station l Barotropic prediction —'— Sx - I Kilometer • 8x = 5000 feet o measure d 8 10 12 14 16 March 30, 1973 18 20 22 24 Houn Figure 33. Comparison of predicted velocities with observations on March 30, 1973 at stations 1 and 2. 102 Computing Centre. The barotropic equations (43) and (44) were evaluated throughout the estuary in a regular advancing computation based on St = 60 seconds with the appropriate boundary conditions. Next equations (45) and (46) were used in turn to update the lower layer depths and velocities by taking the total water depth, h, from the barotropic pass. Finally the upper layer velocities were computed from equation (47). In the lower part of the Main Arm the total water depth, h, is ref.erred to a level model bottom which serves as datum for the stratified equations. The relationship between model and prototype parameters is defined . in Figure 34 and the three conductivity profilers are shown in their March configuration. Model datum is used in the barotropic computation. As shown in Figure 22 there are three flow states of interest: subcritical flow and supercritical outflow and inflow and following each complete integration pass the flow state was defined by evaluating cf from equation (22) and examining the propagation direction of the internal waves at each velocity grid line. If the flow was subcritical everywhere (cj > 0; cJ < 0) the interfacial heights were imposed at each end of the stratified computation domain. The salt water height at the toe was determined by interpolation between the moving salt front (h• fixed at one meter) and the nearest depth section downstream of the domain boundary. When the flows are supercritical out at the mouth (cf < 0; x > 0) the interfacial boundary condition h» at x = 0 is J L z = l6.84m(Geodetic datum) z = l4.24m(Sandheads datum) z = 4.25m (Model bottom) Deepest Channel Bottom z = 0 ( Model datum) J L J 8 10 12 14 16 18 20 22 24 26 Kilometers Figure 34. Relationship between the estuary and model parameters. Conductivity profilers at stations 1, 2 and 3 are shown in their March configuration. o 104 released. Physically this means that the head of salt water can' no longer be maintained at the seaward boundary since flow conditions in the river are determining the water column structure throughout the calculation domain. Somewhere seaward of x = 0 a control section will occur marking the boundary between the river dominated flows and the normal stratified conditions in the Strait of Georgia. Since the characteristic curve corresponding to the positive internal wave no longer enters the finite difference domain from the seaward side (Figure 22) information cannot be drawn from this boundary to advance the solution. Therefore h' at x = 0 is not imposed. Supercritical inflows occured only once in the computations and were confined to the velocity sections near the mouth. Throughout the lower estuary the entire velocity field was flooding upstream and the upper layer flows were supercritical with respect to interfacial waves. Flows were suhcritical in the lower layer and near the toe in the upper layer implying a transition in the upper layer. The unsteady flow conditions would probably produce weak surge waves on the interface propagating upstream from the mouth, and in the numerical scheme it was assumed that the energy lost in such waves could be neglected. The solutions may well be in error for this reason during supercritical inflows, however the overall simulation appears reasonable. The additional internal characteristic entering the domain from the seaward side (Figure 22) requires a velocity boundary condition, supplied by evaluating u using the barotropic flux in eguation (47). 105 Initialization and Wash-Oat The observations imply that the salt wedge can be completely washed out past station 1 and the model must include this feature. Since six grid lines are reguired for an advancing computation wash out is defined, as far as the model is concerned, as occuring when the toe position has retreated downstream of section 6.' The wedge removal on each large ebb tide also means the stratified computation must be re-initialized on subsequent floods. In the more usual tidal models, and by implication in the other stratified flow studies, initialization can be done fairly crudely and if the numerical scheme is convergent, it will overcome the effects of bad starting values and eventually approach the correct solution. Normally these models are "run-in" for two or three tidal cycles to eliminate the initialization effects. In this case the starting values for the interfacial height, h1, are calculated from eguation (28) using a measured value of h' at x = 0. The required values of h are computed by assuming a uniform surface slope between grid lines 2 and 26, u' is set equal to zero and u is calculated from the barotropic discharge. The initialization procedure requires criteria against which the validity of the starting values can be judged. For example, if the surface slope becomes too small the intrusion will be unacceptably large and vice versa. If the initial penetration length exceeds 7 kilometers (grid line 8) or falls below 5 kilometers the starting values are rejected. Further if 106 the value of h» at the wedge toe exceeds two meters initialization fails. By ensuring these restrictions are met, what can only be judged intuitively as reasonable beginning values are obtained. In practice after such a wedge was inserted with u' equal to zero, the lower layer would start to move downstream before reversing to the flood, suggesting that initialization was not too late into the flood cycle. Since this entire procedure was repeated at least once within every 25 hours of computation, the effects of bad initialization are felt in the solutions. The overall validity of the method can only be judged by comparing the predictions to the observations for those parts of the tidal cycle when salt water has penetrated more than 10 or 12 kilometers. It is reasonable to expect the time of arrival at profiler 3 to be affected by bad initialization. The wedge toe has been arbitrarily defined as coinciding with the position where h* =1 meter. In this computational scheme the toe position has been advanced or retreated assuming that the flux of salt water just downstream of the toe, will all go into increasing or decreasing the storage associated with the moving front during each computation interval of 2 St. This assumes that the salt water discharge is constant for 2 St. CoBiparincj Predictions with Observations The basic problem in applying the numerical algorithm centred on finding the correct balance between the two stress 107 terms T± and "Th. Large imbalances affected three aspects of the solutions: the phase and duration of salt water at station 3, the greatest penetration distance during a particular tidal cycle, and the magnitude of the layer velocities. Optimum values for Ki and Kb were found by systematically altering the two coefficients and comparing solutions to data from both measuring periods. The stress formulation in equation (40) gave the best overall fit with final coefficient values of Ki=0.0075 and Kb=0.0055 for the £§rticular barotropic solution discussed previously. Two-layer results The numerical solutions are conveniently compared to observations using the two-layer parameterization of the profiler data and the velocities recorded on March 29 and 30th. In Figure 35 the solutions for February 10 and Harch 17 have been plotted together with the salt layer thickness observations relative to model bottom. The observations Contain errors that become appreciable when the salt water is low on the probe line; this effect is worst at station 3, particularly in the February data. The predicted interface goes to zero almost simultaneously at both stations 1 and 2 due to wash-out as defined earlier. In both solutions the maximum salt water height is reasonably well modelled except for the second flood on March 17th. The phase agreement is satisfactory at station 3 for both records, where the peak intrusion agrees within 40 minutes of the observations. Observations 1 x Station I 0 Station 2 + Stalion 3 MARCH 17, 1973 Hours Figure 35. Comparison of typical interfacial positions predicted by the stratified flow model with two-layer parameterization of observations at stations 1, 2 and 3. 109 In all cases the model tends to flood too early at station 2. The observations near 0000 hours on February 10th and 1700 hours on March 17th show an interfacial separation between stations 1 and 2 averaging about 2 meters, corresponding to a slope of 0.7x10~3. The numerical model could not reproduce such large slopes near the seaward end except .immediately following initialization, when the solution is least accurate. The contour charts in Figures 13 and 14 suggest that mixing of fresh water into the salt layer during reversal of the upper layer velocities is considerable and would have the effect of lowering the interface compared with the no-mixing solution. The model solutions also diverge from the measurements when the ebb tide flow state changes to supercritical at station 1. At this point the boundary condition is released and: the numerical model behaves consistently. The interfacial heights between stations 1 and the vicinity of Canoe Passage junction (x=9 kilometers) fall rapidly until the interface is nearly horizontal and located at a height depending on the penetration length of the toe. Near the toe the interfacial slope increases rapidly forming a rounded nose for supercritical flows along the whole wedge. In the February 10th solution the freshwater layer remains relatively thin (4.1 meters at turn to supercritical as predicted by the model) and the salt water height falls more slowly than the measured response. This is due to the large volume of salt water which must advect out of the estuary in the model. On the other hand with the shorter penetration lenghts during March, the predicted interface falls away very rapidly 110 from the measurements, and provides a poor simulation of the wedge behaviour during this phase. To check this aspect of the Harch predictions, the February 8th solution was computed since it had similar tidal conditions. The same trend occurred but the turn to supercritical outflow lead the observations by only 25 to 30 minutes compared with 50 to 60 on March 17th. The loss in salt water thickness was also more closely matched to measurements. It seems likely the differences are associated with the barotropic model. The spike in the station 1 ebb curve on March 17th resulted from the turn to supercritical outflow conditions. At this point the salt water heights near station 1 collapsed to a lower level and the resulting adjustment of the layer velocites promoted subcritical. flow, allowing imposition of a measured boundary condition some 90 or 100 centimeters above the previous level. Consequently supercritical flows were re-established and remained until wash-out. The velocity data from March 29 and 30th are compared with predicted upper layer velocities in Figure 36. As with the conductivity profilers, velocity measurements were not made in deep enough regions of the estuary to provide a good picture of the lower layer flows. To compute the observed upper layer velocity, u0, the interfacial height at each applicable time was passed to a subroutine which integrated the measurements over the Ill March 29, 1973 — prediction section 5 O observation station 2 March 30,1973 .... prediction section 3 — prediction section 5 A observation station I O observation station 2 8 10 12 14 16 Time ( Hours) 18 20 22 24 Figure 36. Comparison of predicted and measured upper layer velocities at stations 1 and 2 for March 29 and 30, 1973. 112 upper layer; that is, u0 = y7)Ji u(z)dz (49) h where u (z) = measured velocity at height z. It is more difficult to compare the velocities this way than layer depths since the velocity phase varies with distance along the estuary and is computed on a 2 kilometer grid. Station 1 is compared with section 3 and station 2 with section 5 which reflects a separation of only two kilometers in the predictions compared with 2.8 kilometers in the estuary. The measured and predicted velocities agree within ±15 cm/sec for both stations. However, the model solutions appear to lead the observed flows especially at station 2. This is due in part to the distance from the boundary, but is also consistent with the interfacial response. Qualitative Comparisons The numerical solutions for three days in each measuring period have been superimposed on contour charts of constant density in Figures 37 and 38. The densities in terms of °j were derived from the same conductivity data contoured in Figures 13 and 14 using the background temperature data. As before the horizontal axis is time and the vertical columns of zeros indicate the operative probes. The vertical axes give distance above model bottom in meters and show clearly the extent of water column monitored in relation to the assumed river bed. The solid interfacial line at station 1 is the imposed seaward boundary condition and follows exactly the two-layer parameterization except during periods of supercritical outflow which have been o 2 o g u> o O o o o »s ^ CO o . o . o o o CO o 10 o o inittobzation • m«osurad boundary condition STN.I c,'<0 STN.2 Model Bottom FEBRUARY 1973 Figure 37. Interfacial solutions from the stratified flow model are shown superimposed on contour charts of constant density ( Ct) for three days+ in February 1973. The periods of supercritical outflow are indicated by c~ <0 on the station 1 contour chart. r—1 LO o J I MARCH 1973 Figure 38. Interfacial solutions from the stratified flow model are shown superimposed on charts of constant density ( Ct) for three days in March 1973. Periods of supercritical outflow are indicated by <0 on the station 1 contour chart, supercritical inflows by cf >0. M 115 indicated by cj <0 on the charts. The two-layer measurements have also been indicated during supercritical outflow. These figures show that the model is capable of using the barotropic solution to initialize the stratified computation, and with the measured internal boundary condition at the river mouth, calculate the salt water flows for several kilometers upstream. The model provides a reasonable estimate of the duration of salt water at each station and allows confidence in the predicted penetration lengths. As noted previously, the model is least accurate in simulating the layer thicknesses once supercritical outflow occurs. This effect is worst for tides of small diurnal inequality as shown in Figure 38. The numerical model allows the salt water to be washed out of the estuary above Steveston, and judging by the contour charts at station 1 the timing appears accurate within ±30 minutes. The solution for February 10 has also been plotted in a series of longitudinal sections in Figure 39. Each section corresponds in time to the salinity profiles of Figure 12 and by interpolating the salinity data, isohalines have been superimposed on the predictions. The sections also show the computed layer velocites at station 2 and near the wedge toe. Until the turn to supercritical the interface lies near the 15%o isohaline, after which time it appears to rise in relation to the salinity structure. In the numerical model all of the salt water must flow out of the estuary within the lower layer and a further 150 minutes was reguired beyond the final section in Figure 39 116 x ( kilometers ) x ( ki lome te r s ) Figure 39. Interfacial solutions for the tidal cycle on February 10, 1973. Isohalines interpolated from the data in Figure 12 are shown in each section. 117 for this to occur. In contrast, turbulent mixing between the water masses in each layer plays an important role in the estuary, giving rise to a different removal mechanism than occurs in the model. These sections illustrate the behaviour of the numerical model on ebb - the tendency for the interface to level off downstream of station 3 and an increasing slope upstream. I found that the numerical model was influenced by the junction at Canoe Passage which reduced the Main firm flow about 10 per cent on ebb. This effect could account for a change in the interfacial slope due to the decreased velocity and although solutions computed without Canoe Passage exhibit a less pronounced change in slope, they retain the same overall behaviour. This contrasts the limited data in Vreugenhil's paper (1970) which shows the salt layer forming a distinct wedge shape even near maximum ebb although there is no comment on the flow state. One piece of data I was unable to collect was the maximum excursion of the wedge and this would have been valuable. The March tides produced penetration lengths averaging 18 kilometers above Steveston (station 1) while those in February were about 27 kilometers. The very large ineguality on February 11 gave the greatest excursion (36 kilometers) which would put salt water above New Westminster (station 4). Salt water fish species and crabs which are occasionally caught near New Westminster provide indirect evidence to support such a long penetration. I would 118 expect the two-layer model to overestimate the penetration lengths due to the diluting effect of mixing on the salt water layer in the estuary. As far as modelling the gross hydraulics the numerical scheme seems satisfactory, judged ty the duration of salt water and the velocity predictions. It does not perform well on detailed aspects, especially the turn to ebb and maintenance of the salt layer thickness during supercritical outflow. The early flood and ebb characteristic of the solutions, which becomes less severe upstream, seems to be shared by Vreugdenhil's two-layer model (Vreugdenhil (1970), figure 7/) where the predictions appear to peak approximately 40 to 60 minutes too early. Sensitivity Analysis The sensitivity of a number of characteristics of the model solutions for February 10 and March 29 was investigated by holding each stress parameter constant in turn, and varying the other through a range of values. There is a definite limit to values which can be assumed by the interfacial stress Ti since the boundary condition is fixed. If the stress becomes too large, the initialized wedge is repeatedly washed out until the initialization procedure fails. The interfacial stress also provides dissipation of numerical errors and cannot assume a zero value. A mean velocity shear, U=u-u', of 0.300 m/sec is typical of the stratified flows and corresponds to a stress of 2.44 119 dynes/cm2 for Fi=0.36 C7 =h'=6 meters) and the coefficient Ki=0.0075. Numerical solutions were computed for ri+50%, ri + 100% and ri + 30055. A salt water velocity, u'=0.300 m/sec, is also typical and for Kb=0.0055 represents a bottom stress of 4.95 dynes/cm2. Three values were investigated: Tb±50% and Tb=0. The penetration lengths for February 10 are almost independent of Ti, showing a variation of only ±6 per cent about the standard solution (as I shall refer to the solution described in the previous section). On the other hand a zero bottom stress increases the excursion by 40%, and provides a poor simulation at station 3. It is clear from the form of the solutions that Tb is an important parameter. The penetration is more sensitive to changes in Kb than Ki. An important characteristic of the solutions is the time at which supercritical outflow first occurs. In the February 10th solutions this time varied only ±6 minutes for the range of ^"i. However, the zero bottom stress delays the turning point 37 minutes. The effect of this on the layer thicknesses at station 1 is confined to about 5 per cent. Compared with the March solutions, the supercritical outflows on February 10th do not produce such significant errors in the salt layer depth, due possibly to the long penetration, and this may explain the apparent insensitivity to the stress parameters. Both stresses affect the wash-out times at station 1; in 120 this case, the interfacial stress becomes important. With ri+300% the removal time is reduced by 60 minutes, as it is with Tb-5Q%. By comparison the zero bottom stress only delays wash-out 35 minutes compared with about 60 minutes for Ti-503L The removal time and rate appear most sensitive to ^"i but lack of data prevents a conclusive analysis. Some general trends are apparent in the various solutions. Increases in Ki do not produce proportionate changes in TL until the very large ebb currents. Rather the average interfacial stress level remains fairly similar (+10 per cent) but the increased coupling between the layer velocities has the effect of reducing u-u1. Since h' is identically imposed in each solution, the layer depths are not greatly affected, except at station 3. The standard value of Ki (0.0075) seems to lie near a lower limit. Lower values produced marked changes in the wash-out time and the flood-ebb behaviour at station 3, whereas higher values of about the same percentage increase did not. The increased coupling for higher values of Ki has the effect of slowing the upper layer velocities (in fact, collapsing them back onto the barotropic solution) while decreases in Ki allow the upper layer to flood and ebb with higher velocities. The velocity phase is not noticeably changed. Increases in the bottom stress delay the ebb at station 3 and on this basis it is possible to reject values of Kb above 0.007. Variations of ±50% in rb have almost no effect on the upper layer velocity compared with the measurements. As Tb 121 decreases the peak ebb and flood currents in the upper layer are reduced since u' is increased. The reduction is on the order of 7 to 10 per cent for  Tb-50%. The phase appears unaffected. Based on the fitting between observations and the model discussed in the previous section limits can be estimated for the stress coefficients. These values will reflect both the accuracy and extent of field monitoring and the influence of the barotropic model. A minimum value of 0.006 for Ki and a maximum of 0.010 are appropriate. The bottom stress would occupy a smaller range: Kb=0.0055±0.0010. Alternate Interfacial Stress Forms Each of the three alternate stress forms, equations (38), (39), and (41) were substituted into the numerical model and solutions for March 17th were evaluated. The stress coefficient was calculated to give the same stress as equation (40) for u-u^O.300 m/sec and Fi=0.362, equivalent to the values used in the sensitivity analysis. This does not necessarily represent the optimum value for Ki in each form, but was designed to give the solutions a common basis. Superposition of longitudinal sections from each solution revealed that the overall form of the wedge did not differ greatly from the standard. On both flood and early ebb, the largest difference in salt water thickness was 30 centimeters located 10 kilometers above the boundary at station 1. The stress forms providing the weakest coupling, KiyOFi and Kip U | U |, 122 did not allow the wedge to be removed as rapidly as the other forms, and in general were associated with greater lower layer depths. The penetration length predicted by each form varied only ±4 per cent from the standard, but the stresses providing greater coupling delayed the turn to supercritical outflow a total of 56 minutes. The solution for Kip Fi went supercritical first, but the lower stress levels resulted in the longest duration at all stations. Conversely, supercritical flow was delayed the longest for Kip U|U|Fi2 but the removal of salt water was the most rapid once the change of state had taken place. The velocity field is also affected by the degree of coupling; the stress forms in equations (40) and (41) produce the lowest upper layer velocities and greatest salt water currents. The stress form does not change the phase of the velocities, only the relative magnitude of u and u'. During the development of the numerical procedures all stress forms were tried and calibrated. Overall eguation (40) was the most satisfactory, however each stress form permitted a numerical solution which could be fitted to the observations. Equations (40) and (41) are the most responsive to changes in U and with one exception, prevented supercritical inflows. They also produced the greatest delay for the turn to supercritical outflow. I have chosen KyOU|0|Fi since it provided a slightly closer match to observations than the the alternative forms and appears to have reasonable theoretical validity. Unfortunately, 123 the observations coupled with the use of the barotropic solution to provide the upper layer flows and surface slopes (entering the pressure gradient term in the lower layer) prevent a more definitive statement concerning a formulation for the interfacial stress term. Essential Features of The Model A numerical model has been described in this chapter which, although based on some fairly broad assumptions, succeeds in simulating several important characteristics of the observed salt water flows in the Fraser estuary. The model neglects lateral variations in density and velocity, assumes the river can be approximated by a channel of uniform width and depth and draws data required for pressure forces and boundary conditions from a one-dimensional hydraulic model. A further boundary condition is supplied by measurements. Finally, the model neglects mixing of salt and fresh water across the interface, but includes the turbulent stresses using an empirical relation. Consistent with observations, the model predicts that salt water will be flushed out of the estuary on each large daily ebb, even for low flows, and is capable of re-initialization following wash-out. The extent of penetration is predicted within 10 percent of the total depth and the phase shews agreement within ±40 minutes at approximately one-half of the average excursion distance. The model also provides estimates of the average layer velocities which agree within ±15 cm/sec of the 124 observations. On the other hand, the simulation of the layer depths and removal mechanism once supercritical ebb flows are reached is not good. This is to be expected in a model without mixing effects. The shape of the calculated salt wedge is relatively insensitive to both the magnitude and form of the interfacial and bottom stresses. In addition, the phase of the velocity .fieid appears to be almost independent of the stress assumptions. These results are due to the over-riding influence of the barotropic model and to the use of a measured boundary condition. I found during development of both the stratified and barotropic models, that the phase of the layer velocities and depths depend upon the accuracy of the barotropic simulation1 and the estuary geometry. The one-dimensional model described previously represents about the best fit obtainable for the particular schematization. Further improvements could be expected with a better representation of the delta and storage effects along the south bank of the Main Arm. The price for this would be a more complex and expensive computer model. Such improvements, to be properly done, warrant a careful and detailed monitoring programme in the estuary emphasizing velocity measurements. In summary, it is essential to the success of the stratified calculation that the barotropic model be calibrated as accurately as possible. Currents in the upper fresh water layer of the estuary 125 are affected by the salt water flows due to turbulent mixing along the interface. To calculate the correct layer velocities both the relative magnitude and the total value of the interfacial and bottom stresses must by well modelled. A sequence of salt wedge solutions for March 17th is plotted in Figure 40 together with the time variations in both stresses and u-u'. The most important observation is that for nearly all phases of flow the bottom stress, which follows the variation in u*, is larger than the interfacial stress. Also i remains fairly constant during supercritical outflows. This type of variation makes it difficult to generalize, but during ebb the bottom stress is about 60 to 80 percent larger than ^~i. In addition, the maximum unsteady interfacial stresses are considerably greater than those calculated for an arrested wedge. For example, the stationary wedge analyses gave interfacial stresses of approximately 2 to 4 dynes/cm2 compared with values of 10 to 12 dynes/cm2 for the unsteady case. This can be explained in terms of the increased Reynold's stresses near the interface due to turbulent flows in each layer. The interfacial height boundary condition at the river mouth is the forcing mechanism which drives the salt water layer upstream, and the accuracy of the simulation throughout the estuary depends upon the accuracy with which this boundary condition is specified. In two previous studies both Boulot et al. (1967) and Vreugdenhil (1970) have formulated the internal downstream boundary condition in terms of river flow parameters. 126 Distance above Station I ( kilometers) Figure 40. Tidal variations in the interfacial and bottom stresses and the mean layer velocity shear at station 2 for March 17, 1973 (upper). Six longitudinal sections are shown for flood and ebb phases (lower). 127 Basically they postulate, as in the stationary wedge analysis, that a control section will form near the river mouth and that the fresh water layer thickness there will be determined by hydraulic conditions within the estuary. In the presence of an outflowing upper layer c"| = 0 is assumed. Vreugdenhil goes further and assumes that during supercritical inflow both internal waves are stationary and by setting cj" = 0 and solving equation (22) for h' and u' using equation (47): h»(0,t) = h(0,t)/2 + q(0,t)/(2 J g e h (0 , t) ) (50) u'(0,t) = g(0,t)/(2h(0,t)) + / g e h (0,t)/2 (51) where q = the net discharge per unit width. This assumption appears to follow from the work by Schijf and Schb'nfeld (1953) in which they define a "double critical" exchange flow, where the discharge is equal in each layer but opposite in direction. Then for unique flow conditions the upper layer velocity tends to arrest one internal wave and the lower layer flow arrests the other, producing a case where ct -c~—-*• 0. Rigter (1970) has suggested that such a situation gives the maximum discharge in the salt water layer corresponding to optimum hydraulic conditions. The generality of this kind of assumption must be questioned since during calculations for the Fraser River, the tendency to supercritical inflow corresponded with upstream velocities in each layer and no mechanism existed for arresting the positive internal wave. These boundary conditions were evaluated for the March 17th tide and compared with the measured values for h1 and the lower layer velocity calculated directly from the momentum eguation in Figure (41). It must be emphasized 128 .5r-f \> Supercritical inflow 2 N boundary condition o ®0.5 E \ v / v / \ \ / \ ^ /\ w x Computed boundary velocity «L \ ./ I 1 1 1 I *^-<S I i JL > 4 6 8 10 12 14 16 18 20 22 24 March 17, 1973 Hours Figure 41. Comparison of the measured boundary condition for h' at station 1 with the theoretical condition of Vreugdenhil (1970) (upper). Comparison of u' calculated from the stratified flow model and the theoretical relation of Vreugdenhil (lower) . 129 that the Fraser model flows did not flood supercritically and Vreugdenhil applies equation (51) only to such a case. The phase agreement between h* from equation (50) and the observed h1 is particularly poor as is the duration of high salt water, however the maximum heights agree within one meter. Equation (50) predicts the maximum h' at a time corresponding very closely to the maximum in g(0,t). Since the phase of both u' from equation (51) and calculated in the model depend largely on the barotropic calculation they are closely matched, but the velocities are consistently 40 to 80 cm/sec apart. Such large salt water flows, used as a boundary condition in the numerical model destroyed the stratified calculation during a series of tests where supercritical inflow was induced by lowering Ti. The doubly critical flow assumption cannot be used to provide two internal boundary conditions for the Fraser estuary since under flood conditions the entire water column flows upstream. It is reasonable to expect that other shallow estuaries with large tides will behave similarly and an accurate specification of flow conditions at the river mouth can be made with the relatively straightforward measuring techniques developed as part of this work. 130 Chapter 6. SUMMARY AND CONCLUSIONS The work carried out and reported in this thesis falls into two distinct but complementary catagories, firstly a programme of field measurements to reveal the nature of the salt water intrusion and secondly a theoretical study into the relationship between tidal and discharge conditions and some important characteristics of the salt wedge motions. The mathematical models also provide a tool for predicting the outcome of changes to the estuary, an aspect which is only touched upon by way of two examples. Two series of observations were carried out in the winter of 1973 for very low fresh water discharges and tides of large and small diurnal ineguality. Self-contained instruments were installed at three points in the estuary which monitored the conductivity structure for periods of several days. Significant salt water intrusions were detected, the extent and duration of which were determined by the nature of the Strait of Georgia tides. Penetrations of 15 kilometers or more were produced by tides of large diurnal inequality having a long high water duration. The salt wedge was also removed from the estuary (upstream of Steveston, B.C.) on each large daily ebb, in spite of the low winter discharges. The conductivity profiles indicate that both salt and fresh water were mixed across the interfacial region and the absence of large density gradients suggested that turbulent mixing took place within each layer. Surface salinities usually 131 exceeding 10 %B accompanied high tide near station 1 and longitudinal salinity gradients were present both near the surface and at depth. These characteristics place the lower Fraser River into the "moderately stratified" class of estuaries (Cameron and Pritchard (1966)), although somewhat of a unique example since it appears to be flushed of salt each tidal cycle. In these estuaries the distribution of salinity is governed by the horizontal advection of salt water along with vertical advection and turbulent diffusion between the layers. The unsteady dynamics of the Fraser estuary have been simulated with a two-layer mathematical model in which the layers were defined in terms of the density structure. The mixing of salt and fresh water across the interface was neglected but the turbulent momentum fluxes were empirically included. The model was solved in explicit finite difference form in conjunction with a barotropic Calculation over the entire estuary. The use of a measured boundary condition for the salt water depth enabled reasonable predictions of the extent and duration of salt water in the lower estuary as a function of the tide at Steveston. The modifying effects of the salt water flows on the fresh water currents were predicted by the model. The success of the model in predicting the phase and magnitude of currents and penetration depended primarily on the accuracy of the barotropic calculation. This was also Vreugdenhil's conclusion for a similar model of the Rotterdam Waterway (Vreugdenhil (1970)). Of nearly equal importance in 132 terms of the duration of salt water throughout the estuary was the measured boundary condition for the interface at the river mouth. The duration of salt water and current magnitudes were also a function of the correct balance between the bottom and interfacial turbulent stresses; typically the interfacial stress had maximum values on the order of 10 dynes/cm2, and the bottom stress was usually about 80 percent greater, namely 18 to 20 dynes/cm2. Both stresses are significant in the dissipation of energy in the flows. Four empirical relations for the interfacial stress were substituted in the numerical model, expressing increasing degrees of sensitivity to the mean velocity shear between the layers. Due to the rigid boundary condition for the salt water depth and the resolution of the data, it was not possible to make a conclusive choice among them. Each relation could be calibrated to the data in much the same form and the final choice of KioU|U|Fi was made mainly on the correspondence to the measured velocities. The interfacial Froude number Fi has been included in the stress formulation to represent the influence of increased mixing near the onset of supercritical flows. The bottom stress was modelled using the quadratic relation KbyOu'|u'| and the calibrated values for each coefficient were Ki=0.0075 and Kb=0.0055. The predicted interfacial positions were relatively insensitive to small changes in these coefficients; differences were seen in the magnitude of the layer velocities. Using mainly 133 the correspondence between measured and predicted fresh water velocities and the interfacial predictions at station 3 the following ranges for each coefficient were estimated: 0.010 < Ki < 0.006 and 0.0065 < Kb < 0.0045. If the stress coefficients were within these ranges the model predicted salt water thicknesses within 10 percent of the total water depth and a phase agreement at maximum intrusion within ±40 minutes compared with measured values obtained approximately 10 kilometers upstream of the mouth. The fresh water velocities showed agreement within ±15 cm/sec of the observations. Maximum ebb currents were on the order of 2 meters/sec. £2l!£ludin,g Remarks Two theoretical studies of unsteady salt wedges have been reported in the past, both of which are based on two layers and neglect the mixing of salt and fresh water across the interface. In only one case were the predictions compared with observations and these measurements occupied less than one tidal cycle. Both models were solved in explicit difference form using Lax-Wendroff techniques. The field monitoring programme forming part of the present work succeeded in producing a detailed description of the salt wedge motions in the Fraser River over several days and different tidal conditions. This allowed a rigorous comparison to be made between the predictions from a two-layer "no-mixing" model and the actual hydraulic behaviour in the estuary. I found 134 such a model provides reasonable estimates of the extent and phase of the penetration, as well as the currents in each layer if the correct boundary conditions are applied. In this study the interfacial height at the river mouth was found from measurements. The model does not provide a good simulation of the layer thicknesses near the river mouth during supercritical outflows, however even during this flow state the upper layer velocities are well modelled. Two examples are discussed in Appendix A which illustrate some of the more important aspects of the salinity intrusion as it relates to sedimentation and water quality. Basically the hydraulic models serve to provide information on the variations in the salt water location and the velocities in each layer throughout the estuary, information which is reguired to make guantitative predictions for sedimentation or water guality through the use of dispersion models. Some of the more important aspects of the salt water intrusion found in the examples are: (i) Dredging of the shipping channel to provide an additional 3 meters of clearance would increase maximum penetrations about 20 percent. This would result in salt water located in and above the trifurcation at New Westminster for much greater periods of time than at present and would increase the sedimentation in this region which is already dredged on a continual basis. (ii) The time required to flush pollutants confined to the fresh water layer is significantly decreased by the salt wedge presence compared with times computed using the unstratified velocity field. (iii) The distribution of parcels of contaminated water within the estuary is significantly changed by the baroclinic flows compared with the unstratified distribution. In general, the stratified computation shows that water parcels originating both from Gilbert Road and Annacis Island pass less often through outfall regions, which would tend to increase effluent loading on 135 the fresh water layer, than those water parcels computed with the barotropic velocities. (iv) Host high tides during the measuring periods produce predicted excursions exceeding 18 kilometers above station 1, thereby bringing stratified flows into the vicinity of the Annacis Island sewage treatment plant. The suppression of vertical mixing associated with the stratification means that effluents discharged into either layer might tend to remain in that layer resulting in higher concentrations than would be the case in unstratified rivers. -(v) Increased salt water penetration associated with further dredging would increase the duration of stratified flews around Annacis Island, thereby increasing effluent concentrations due to reduced rates of vertical mixing. It is clear that the salt wedge has important effects on circulation in the estuary and that subsequent hydraulic models which may be used for predicting changes in sedimentation or water quality, must include the baroclinic flows. Recommendations for Future Research The numerical model for stratified flows discussed in this thesis neglects the effects of salt and fresh water mixed across the interface. The effect of longitudinal density gradients on the equations of motion is usually negligible (Rossiter and Lennon (1960) ) however the variations in the density contrast along the wedge can be expected to influence the interfacial stress when parameterized in terms of a Froude number. Therefore the logical next step in the development of the hydraulic model is to incorporate mixing across the interface, formulated in terms of an entrainment term plus a turbulent diffusion term. The main difficulty involved in this step is the introduction of two additional empirical coefficients 136 for the mixing terms, a difficulty which can be resolved however, by obtaining -sufficient data from the estuary to derive the coefficients from a salt balance calculation. Mixing effects were not included in the present study since they introduce considerable complexity to the numerical model and make extensive demands for field data, both to determine the mixing coefficients and to properly verify the predictions. For example, the observations on March 30, 1973 were insufficient to permit determination of the mixing rates, due in part to the loss of probes at stations 1 and 2 and in part to the location of each station. In rivers like the Fraser, a proper salt balance calculation will probably require enough stations across the width to provide a cross-sectional average of the currents and salinities due to secondary currents produced by the bends. Three stations across the width, one located near the maximum depth and the other two located mid-way between the first station and the banks would probably be adeguate, but one can immediately appreciate the problem of logistics in such a project! Six stations operating simultaneously with at least four located in the shipping channel. There is the additional consideration that a mixing model also requires the specification of the salinity or density at the river mouth, although this can be obtained by a suitable parameterization of the type of data collected in this project. Calibration of numerical models is necessary because we need to include empirical relations for turbulent stresses or diffusion terms and in this regard, velocity data from the 137 prototype is particularly valuable. Furthermore this data should be collected across the width of the estuary where feasible to allow cross-sectional averaging. If only one station can be occupied it should be located as close as possible to the deepest part of the channel. A number of recommendations can be made for future monitoring programmes similar to the one reported here. The problem of conductivity probe loss could be significantly reduced by casting the probe coils in plastic or another waterproof compound rather than the hollow sealed units used in this project. The three pile triangular platforms were satisfactory as designed but could be made more durable by increasing the pile clamp thickness from 1/4 inch to 3/8 or 1/2 inch. The clamps at station 1 failed through fatigue along the weld between the clamps and the steel tube fame. To be certain of eliminating or maintaining constant wire angles on the instrumentation, submerged anchor weights exceeding 300 pounds on the guide wires are recommended. 138 BIBLIOGRAPHY Abbott, M., 1961. "On the Spreading of One Fluid Over Another, Part I", La Houille Blanche, Vol. 16, No. 5, p. 622 - 628. Abbott, M., 1961. "On the Spreading of One Fluid Over Another, Part II", La Houille Blanche, Vol. 16, No. 6, p. 827 - 846. Abbott, M., 1966. An Introduction to The Method of £!}2l§£i2£isticsx Thames and Hudson, London. Baines, W., 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, Ottawa. Baines, W., 1953. "Survey of Tidal Effects on Flew In the Fraser River Estuary, June 10 and 11, 1952", Report No. MH - 40, National Research Council of Canada, Ottawa. Boulot, F., Braconnot, P. and Marvaud, Ph., 1967. "Determination Num^rigue des Mouvements d'un Coin Sale", La Houille Blanche, Vol. 22, No. 8, p. 871 - 877. Cameron, W. and Pritchard, D. (1963). "Estuaries", Chapter 15 in The Sea (ed. M. Hill), Vol. 2, John Wiley & Sons, New York, p. 306 - 324. Courant, R., Friedrichs, K. and Lewy, H., 1928. "Uber die Partiellen Differenzengleichurgen der Mathematischen Physik", Mathematische Annalen, Vol. 100, p. 32 - 74. Dronkers, J., 1964. Tidal Computations in Rivers and Coastal Waters^ North Holland Publishing Company, Amsterdam, The Netherlands. Dronkers, J., 1969. "Tidal Computations for Rivers, Coastal Areas, and Seas", Journal of the Hydraulics Division, ASCE, Paper No. 6341, HY1, p. 29 - 77. Ellison, T. and Turner, J., 1959. "Turbulent Entrainment in Stratified Flows", Journal of Fluid Mechanics, Vol. 6, p. 423 - 448. Farmer, D. and Osborn, T., 1973. "An Instrument for Measuring Conductivity Profiles in Inlets", Journal of the Fisheries Research Board, Canada, 29(12), p. 1767 - 1769. Farmer, H. and Morgan, G., 1953. "The Salt Wedge", Proceedings of the Third Coastal Engineering Conference, p. 54 - 64. Grubert, H. and Abbott, M., 1972. "Numerical Computation of Stratified Nearly Horizontal Flows", Journal of the Hydraulics Division, ASCE, Paper No. 9300, HY10, p. 1847 - 1865. 139 Hodgins, D. and Quick, M., 1972. "Computer Studies of Estuary Water Quality", Proceedings 13th Coastal Engineering Conference, July 10 - 14, 1972, Vancouver, Canada, p. 2327 - 2338. Joy, C. , 1974. Water finality Monitoring in Estuaries^, Ph.D. Thesis, University of British Columbia, Vancouver, Canada. Keulegan, G., 1966. Chapter 11 in Estuary and Coastline Hydrodynamics (ed. A. Ippen), McGraw-Hill, New York, p. 546 -574? Lafond, E., 1951. EE2£§§§iS9. 22§i.nographic Catax Hydrographic Office Publication No. 614, U. S Navy Hydrographic Office, Washington, D. C. , p. 9 1. Lofquist, K., 1960. "Flow and Stress Near an Interface Eetween Stratified Fluids", The Physics of Fluids, Vol. 3, No. 2., p. 158 - 175. Miles, J. , 1961. "On The Stability of Heterogeneous Shear Flows", Journal of Fluid Mechanics, Vol. 10, p. 496 - 508. Pritchard, D., 1956. "The Dynamic Structure of a Coastal Plain Estuary", Journal of Marine Research, Vol. 15, No. 1, p. 33 - 42. Rattray, M., 1964. "time-Dependent Motion in an Ocean; A Unified Two-Layer, Beta-Plane Approximation", Studies in Oceanography^ (ed. K. Yoshida), University of Washington Press, Seattle, Washington, p. 19 - 29. Richtmyer, R. and Morton, K. 1966. Difference Methods for Initia1^value Problems^ 2nd edition, Interscience Publishers, New York7 p.~288. Rigter, B., 1970. "Density Induced Return Currents in Outlet Channels", Journal of the Hydaulics Division, ASCE, Paper No. 7086, HY2, p. 529 - 546. Roach, p., 1972. Computational Fluid Dynamics^ Hermosa Publishers, Albuguergue, New Mexico, U.S.A. Rossiter, J. and Lennon, G., 1965. "Computation of the Tidal Conditions in the Thames Estuary by the Initial Value Method", Proceedings of the Institution of Civil Engineers, Paper No. 6855, p. 25, 56. Schijf, J. and Schonfeld, J., 1953. "Theoretical Considerations on the Motion of Salt and Fresh Water", Proceedings Minnesota International Hydraulics Convention, I.A.H.R., September 1-4, 1953, Minneapolis, Minnesota. Smith, G., 1965. Numerical Solution of Partial Differential J. Oxford University Press, Oxford. 140 Stommel, H. and Farmer, H., 1952. "Abrupt Change in Width in Two-Layer Open Channel Flow", Journal of Marine Research, Volume XI, No. 2, p. 205 - 214. Taylor, G., 1931. "Internal Waves and Turbulence in a Fluid of Variable Density", Conseil. Perm. Intern. Pour L'Expl. De la Mer, Rapp. Et Proc-Verb., Vol. 76, p. 33 - 42. Turner, J., 1973. Buoyancy Effects in Fluids^ Cambridge University Press, Cambridge. Vreugdenhil, C, 1970. "Two - Layer Model of Stratified Flow in an Estuary", La Houille Blanche, Vol. 25, No. 1, p. 35 - 40. Waldichuk, M., Markert, J. and Meikle, J., 1968. "Fraser River Estuary, Burrard Inlet, Howe Sound and Malaspina Strait Physical and Chemical oceanographic Data, 1957 - 1966", Manuscript Report Series No. 939, Fisheries Research Board of Canada, Pacific Oceanographic Group, Nanaimo, B. C. 141 APPENDIX A. ENGINEERING ASPECTS OF THE STUDY The numerical models developed in this thesis can be used -to provide part of the information required to predict the outcome of changes in the estuary resulting from modifications in the present forms of exploitation. The salinity intrusion produces effects which must be considered in problems dealing with sedimentation and water quality in the lower estuary. Generally speaking the hydraulic models provide the time varying velocity fields in each layer and the extent of salt water penetration. Submodels which relate sedimentation or the dispersion of possible contaminants to the velocity information must be derived and used in conjunction with the basic hydraulic models before a predictive capacity can be realized. However it is possible to illustrate some of the salt wedge effects by considering only the information provided by the hydraulic models. Sedimentation is accentuated by the salt water since it produces a region of turbulent mixing near the toe, slowing down the velocities next the bottom and promoting deposition. This may be accompanied by flocculation which further increases sedimentation. The variations in toe position and salt water duration upstream of the mouth are important parameters. It is of interest to examine the effects of two imagined changes in the estuary: a further dredging of the shipping channel to add three meters to the average depth or, on the other hand, allowing the estuary to return to its natural state (reducing the depth by 142 about 2 meters on average). I have evaluated the model for February 10th using the measured boundary conditions but altering the total depths. The results are summarized in Table 4 in terms of the penetration lengths for three times in the tidal cycle, and the wash-out time. All runs were initialized at the same time. As expected the deepest channel admitted the longest wedge; the shallowest showed the shortest penetration. Since the measured boundary condition was used for all runs, the fresh water depth along various wedges remained similar and the differences in penetration compared with the present mean depth arise very much as an extrapolation of the "standard" interfacial solution. On the basis of these values we might expect additional dredging of two or three meters to increase the penetration about 20 percent over the present excursion distances. A reduction in penetration of approximately 15 percent would accompany depths decreased on the order of two meters. The salt water duration predicted by the model is not significantly altered for changes of this magnitude and this would also be the case in the estuary since the shipping channel usually occupies less than 15 percent of the width. The significance does lie in the changes brought about by deepening. The increased penetration means that over many tidal cycles, the presence of salt water in the trifurcation area at New Westminster will be increased and sedimentation enhanced in this area. This part of the river, which presently receives 143 Table 4 Salt wedge intrusion characteristics in the Fraser River for various model depths. The deep channel is 3 meters below the present and the shallow is 2 meters above the present. ESTUARI TYPE Deep PENETRATION (KILOMETERS) Present Shallow Time Feb. 10, 1973 0300 hours 21.6 17.4 14.3 0500 21.0 16.2 12.9 1200 32.0 27.2 23. 5 Wash-out 17h53m 17h43m 17h27m 144 continual dredging, would be subjected to increased rates of deposition and the revenues generated by expanded shipping would have to be balanced against the cost of channel maintenance. Hater quality is rapidly becoming an issue on the Fraser River and has been investigated recently by both field monitoring and the development of one-dimensional dispersion models (Joy (1974)). In these dispersion studies, baroclinic effects are usually neglected and the validity of this assumption can be examined with the two-layer model, at least as far as flushing in the lower estuary is concerned. Neglecting longitudinal and vertical mixing, let us follow two particles advected by the fresh water flows: one introduced into the two-layer model and one into the barotropic velocity field for a period of two days. Two separate runs have been made using a total of four particles. Numbers 1 and 2 were injected offshore of the Gilbert Road sewage treatment plant on a flooding tide and particles 3 and 4 released near the proposed Annacis Island sewage treatment plant during an ebb tide. The course followed by each particle is sketched in Figure 42 in relation to the Main Arm of the Fraser River. I have assumed the particles will all remain in the Main Arm flows. The main effect of the salt water layer is to increase the fresh water ebb velocities compared with the barotropic velocity field, and this can be seen in the particle 1 and 2 paths. The "stratified" particle (1) advects slightly further upstream on flood but significantly further downstream on the 145 Figure 42. Advection paths of four particles released into the Main Arm of the Fraser River. The small numbers along each path indicate the time in hours following release. 146 subsequent ebb. Also its return is modified and of consequence,' it does not pass back into, the outfall region. Particle 2 does flood a second time into the outfall area and flushing is delayed about 3 hours compared with number 1. Much the same trends appear in the Annacis Island particle paths and it is to be noted that the barotropic model does not flush it s particle until one tidal cycle beyond that of number 3. In general the time required to flush particles in the fresh water is reduced by the salt wedge, the extent of reduction depending upon the point of release. These particle paths imply that there would be considerable differences in predicted concentrations for contaminants between dispersion models based in one case on the stratified velocity field and in the other case, based strictly on the barotropic velocity field. For example, the barotopic particle released at Annacis Island (4) passes five times through the Gilbert Road outfall region compared with three passes for the stratified particle (3). Thus predicted effluent loading on the fresh water near the river mouth would be much higher from a barotropic dispersion model than one including the salt water effects. The same trend is also apparent in particles (1) and (2), both released from Gilbert Road. These results show that the baroclinic flows do have a significant influence on the purely advective part of dispersion and must be taken into account in numerical dispersion models. 147 The Effects of Mixing and Stratification The observations show that turbulent diffusion across the interfacial region occurs along the salt wedge and produces a two-way transfer of water properties. This means that deleterious substances released into either layer will, in some proportion, end up in the other water mass before being flushed from the estuary. A complete simulation of effluent loading or the dispersion of contaminants would therefore require the incorporation of diffusion terms into the dispersion models. The stratified flow model could be used to provide the basic hydraulic information, however, it is clear that once the rates of vertical mixing between the layers have been measured for use in the dispersion studies, the hydraulic model itself can be improved by allowing the layer densities to vary. The salt wedge intrusion is also important from the standpoint of vertical mixing within each layer. In the unstratified river, discharges at either the surface or bottom are usually mixed over the whole depth within one to two hours of release and as a result they suffer considerable dilution. When the river is stratified the rate of vertical mixing is suppressed and effluents tend to be confined to thinner layers with higher concentrations than the corresponding unstratified situation. This is compensated for in some degree by the decreased time required to flush the effluent from the estuary. Predicting these influences on the dispersion of contaminants relies on a knowledge of the salt wedge motions. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items