CHARACTERISTICS OF TIDALLY-FORCEDPOLLUTANT TRANSPORT INNARROW CHANNELSbyMichael J. B. ColeB.Sc.E., Queen's University, Kingston, 1989A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESDepartment of Civil EngineeringWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIASeptember 1992© Michael J. B. Cole, 1992DateDepartment ofThe University of British ColumbiaVancouver, CanadaIn presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)DE-6 (2/88)AbstractPhysical and numerical model studies have been conducted to lend insight to thesubject of transport and dispersion of a neutrally buoyant effluent releasedinstantaneously from a point source into a one-dimensional tidal body. Closed-form expressions are obtained for the fluid velocity and free surface elevationdue to tidal forcing in a one-dimensional channel. Physical experiments wereconducted in the Hydraulics Laboratory of the Department of Civil Engineeringat the University of British Columbia, and measurements of both effluentconcentration and velocity along the test channel were recorded. A numericalmodel based on the Finite Difference Method was developed for computing thetime-varying effluent concentration within the channel. The influence ofgrid size and time step used in the numerical model was examined.The following parameters were varied: effluent discharge location, release timewithin the tidal cycle, average water depth, tidal amplitude, and tidal period. Theexperimental measurements and numerical predictions were correlated to derivean effective one-dimensional diffusion coefficient for the numerical model. Thenumerical results suggest that the fluid velocity at the instant and location ofpollutant release is of primary importance to the mixing of the pollutant cloudinto the surrounding waters.iiTable Of ContentsPageAbstract iiList of Tables vList of Figures viAbbreviations ixAcknowledgment 1 Introduction 11.1 General 11.2 Assumptions and Direction of Research 31.3 Review of Pollutant Transport 62 Theory 112.1 Governing Equations 122.2 Initial and Boundary Conditions 142.3 Development of Hydrodynamic Component 162.4 Dimensional Analysis 193 Description of Numerical Modelling 203.1 Hydrodynamic Component 213.2 Advection/Diffusion Component 213.2.1 Advec tion 253.2.2 Diffusion 293.2.3 Advection/Diffusion 323.3 Combined Components and Extensions 35iii4 Description of Physical Experiments 424.1 Testing Procedure 424.2 Development of Apparatus 434.2.1 Effluent 484.2.2 Aspects of Flow 495 Results and Discussion 515.1 General Observations 515.1.1 Influence of Effluent Discharge Location, x p/X. 515.1.2 Influence of Depth, da 525.1.3 Influence of Tidal Amplitude, A/d 525.1.4 Influence of Tidal Period, T 535.1.5 Influence of Discharge Time, to/T 535.2 Combined Results 535.3 Discussion 586 Conclusions and Recommendations 606.1 Conclusions 606.2 Recommendations for Further Study 62References 64Tables 66Figures.. 72ivList of Tables3.1 Numerical model testing programme for a pure advection case. 673.2 Numerical model testing programme for a pure diffusion case. 683.3 Numerical model testing programme for the advectiveldiffusive componentwith constant velocity. 694.1 Physical testing programme showing variations in tidal conditions based uponXp , T, A, and t/T. 705.1 Best Fit Analysis comparing physical concentration traces with numericallyproduced traces for various diffusion coefficients. 71List of Figures1.1 Sketches of (a) 1-D, (b) 1.5-D, and (c) 2-D channels. 732.1 Definition sketches of one-dimensional situation (a) tidal notation,(b) elevation, and (c) plan. 743.1 Layout of the Finite Difference scheme showing boundary conditionsand numerical cells 753.2 Initial shapes of pollutant concentrations tested in numerical models 763.3 Effect of varying the Courant number (=UAtlAx) within the advectionmodel 773.4 Effect of varying At & Ax while holding the Courant number (=UAtl Ax)constant within the advection model 783.5 Effect of varying y (= DAtlAx 2) within the diffusion model byvarying D while holding At & Ax constant 793.6 Stability range for FTCS Finite Difference operator (after Leonard, 1979) 803.7 Effect of varying Pe (= Cr I 0 within the advection/diffusion modeland varying Cr, At , and Ax while holding 7, U, and D constant 813.8 Effect of varying y(= DAtl Ax2) within the advection/diffusion modelwhile holding At & Ax constant 823.9 Velocity variations along channel length during various stages.ofthe tidal cycle as predicted by closed-form solutions for the Base Case 833.10 3-D perspective plot showing path of pollutant cloud in both timeand space 84vi3.11 Concentration contours in time and space for the Base Case 854.1 Plan and elevation of channel used in the physical model studies 864.2 Set-up of physical model instrumentation 874.3 Position of conductivity probes, x c, in relation to various effluent inputpositions, xp 884.4 Correlation for the conductivity probe between measured conductivity andthe predetermined salinity values 894.5 Video images showing the movement of the pollutant cloud for theBase Case 905.1 Repeatability of effluent concentration measurements from physicaltests for the Base Case 915.2 Concentration traces at locations xc = 3.0, 3.5, and 4.0 m comparingphysical and numerical results corresponding to Test Series 2a 925.3 Concentration traces at locations xc = 3.0, 3.5, and 4.0 m comparingphysical and numerical results corresponding to Test Series 2b 935.4 Concentration traces at locations xc = 3.0, 3.5, and 4.0 m comparingphysical and numerical results corresponding to Test Series 3a 945.5 Concentration traces at locations x c = 3.0, 35, and 4.0 m comparingphysical and numerical results corresponding to Test Series 4a 95vii5.6 Concentration traces at locations x c = 3.0, 35, and 4.0 m comparingphysical and numerical results corresponding to Test Series 5a 965.7 Concentration traces at locations xc = 3.0, 35, and 4.0 m comparingphysical and numerical results corresponding to Test Series 5b 975.8 Velocity profiles comparing physical and numerival results at variouslocations for the Base Case 985.9 Numerical concentration contours comparing various tidal amplitudeto depth ratios, Aid = (a) 0.05, (b) 0.09, and (c) 0.19 995.10 Numerical concentration contours comparing various input positionsalong channel length, xpI A = (a) 0.20, (b) 0.47, and (c) 0.74 1005.11 Numerical concentration contours comparing various input timeswithin tidal cycle, tpiT = (a) ebb tide and (b) flood tide 101ViiiAbbreviationsA = tidal amplitude (m)Ar = cross-sectional area in the y and z directions (m2)c = shallow wave celerity (m/s)Co = initial concentration (mg/1 or parts per thousand )C = pollutant concentration (mg/1 or parts per thousand )Cr = Courant numberd = average water depth (m)D = diffusion coefficient (m2/s)g = gravitational constant (m/s2)H = total water depth (m) equal to d +11k = 2n/L (m-i)K = roughness parameterKd = decay coefficientM = Mass of pollutant (kg)Pe = Peclet numberq = source/sink termRe = Reynolds numbert = time (s)to = initial discharge time (s)T = tidal period (s)Te = duration of effluent outflow (s)U = velocity in the x-direction (m/s)V = velocity in the y-direction (m/s)W = velocity in the z-direction (m/s)x = position along channel direction (m)xc = position of pollutant concentration measurement (m)xP = position of pollutant discharge (m)Y = position in the cross channel direction (m)z = position in the vertical direction (m)• ratio of advection to diffusion11• water elevation relative to the still water level (m)• length of channel in the x-direction (m)fluid density (kg1m3)• = 27JT (radls)• = kinematic viscosity (m2/s)8 Dirac functionAcknowledgmentThe author would like to express his profound gratitude to his supervisor Dr. M. Isaacson forhis help and encouragement throughout this project. Thanks are also extended to Dr. G. Lawrencewho assisted in the supervision of this thesis and offered valuable insight into its subject.The author would also like to acknowledge his friends and colleagues in the Coastal Engineeringgroup at UBC without whom this thesis would not have materialized. The never-ending assistanceof S. Prasad and N. Yonemitsu is especially appreciated. Technical support from Kurt, and Ronof the Civil Engineering technical staff is also gratefully noted.Finally, the author wishes to thank Benjamin B. for providing a most conducive atmosphere inwhich to write the majority of this thesis.x1Chapter 1 - Introduction1.1 GeneralThe characteristics of pollutant transport within channels and harbours are of interest for variousreasons. With increasing urbanization and use of recreational watercraft along the BritishColumbia coast, the amount of pollutants entering our coastal waters is expected to increase.Coincident with this increase is an ever-growing concern for the environment and our waterquality. These concerns affect decisions relating to the design and placement of man-made coastalfeatures such as harbours and outflow pipes and influence the extent and manner in which effluentsare discharged.Local pollutants that enter enclosed waters as a direct result of boating activities include anti-fouling paint, fuel spills, and runoff from various related works (such as loading docks andshipyards). In addition, non-related pollutant sources include industrial effluent discharged fromtributaries, sewage and storm water outfalls. These pollutants affect a channel's water quality andultimately affect the water quality outside the channel. The extent of any environmentalconsequences related to a specific pollutant release depend on:1) toxicity,2) initial pollutant concentrations,3) the rate at which the pollutant is mixed into the surrounding waters, and4) the subsequent rate of decay or rate of expulsion (or flushing) of these affected watersfrom the local zone.The rate of mixing depends on various hydrodynamic factors (including circulation) of thechannel and on the diffusion characteristics of the pollutant. One would expect that the faster the2mixing and flushing occur, the less severe the immediate environmental effects (being related tolower overall pollutant concentrations within the channel).Generalized CaseIn the general case, some of the more significant factors affecting pollutant transport include thefollowing:Pollutant Outflow - point source ( continuous or instantaneous )- line source ( continuous or instantaneous )- multiple sources at various locations within a basin or channel- stage of a tidal cycle at which pollutants are releasedStratification - estuarine conditions of prominent fresh and salt water layers- variations in density between the pollutant and the surroundingwaterTidal Range - variations of the tidal range to average depth ratioReturn Flow - variations in the tidal prism describing the return of pollutant to thechannel on the flood tideHydraulic Forcing - river and groundwater sources causing mixing and/or stratificationChannel Geometries - broad or narrow entrance conditions- length to breadth variations- length to depth variationsExternal Mixing - wind forcing- wave actionHydrodynamics - linear and nonlinear considerationsFriction Effects - variations in flow and turbulenceFlow Constraints - retarding effects of irregular channels and mooring structuresChapter 1 - Introduction3The tidal prism mentioned above, is defined as the percentage of harbour water that exitted on theebb tide and returns from the sea during a flood tide. This percentage is dependent on the effects oflogshore drift.Local Concerns The significant tidal range in British Columbia suggests the occurrence of large flushing andtransport rates in coastal channels and harbours. In channels such as False Creek, Vancouver,where there are complex pier interactions and the channel shape is long and narrow, one wouldexpect flushing rates to be low and mixing rates to be unpredictable. Because of the rockycoastline in British Columbia, groundwater intrusion is usually not a concern. Local winds alongthe coast are irregular and seasonally dependent. Phenomena such as outflow winds can createextremely strong surface shears for extended periods during winter months. Furthermore, sincethe summer months are the period of greatest channel activity (due to pleasure craft), flushing andtransport mechanisms dominant during this period are of the greatest importance for any boatrelated pollutants.Potential Benefits of ResearchIf the transport characteristics of a pollutant are known, monitoring programs can be conductedmore efficiently, by testing at locations thought to be most affected. In the case of pollutants thatneed to be traced at levels close to their detection limits, time and money can be saved by omittingdilute areas from monitoring programs Remediation and worst-case scenarios can be forecast toaid in possible clean-up operations.1.2 Assumptions and Direction of ResearchThe term 'pollutant transport' is used to denote the movement of pollutant due to advective anddiffusive processes. The advective processes move both the pollutant and the surrounding fluidalong the channel without diluting the pollutant or distorting it's shape. The fluid's movement isChapter 1 - Introduction4governed by hydrodynamic processes relating to the effects of viscosity, density, and gravity.The diffusive processes at work within this system are comprised of two main components:1) shear flow dispersion2) turbulent diffusionShear flow dispersion is, itself, a function of the hydrodynamic conditions within the channel.Vertical velocity variations cause shearing of the fluid and molecular diffusion occurs over thisinterface. Turbulent diffusion is a function of the channel velocity and of the turbulent mixinglength which decribes the magnitude of the turbulent processes at hand. The major source ofmixing within the experiment is due to shear dispersion which is typically the major source ofdiffusion in estuaries (Fischer et al., 1979).Pollutant transport in an inlet can be characterized as one, two, or three dimensional on the basisof the importance of depth-wise and width-wise variations of flow velocities and resulting pollutanttransport. A one-dimensional (denoted 1-D) situation corresponds to no significant depth or widthvariations and corresponding low velocities in these directions.In this context, Figure 1.1 shows three separate shallow water channel configurations denoted 1-D, 1.5-D, and 2-D. These three situations are based on the assumption that the channel depth isrelatively shallow so that mixing occurs quickly in the vertical and velocities in the vertical directionare negligible.In the first case (1-D) , it is assumed that mixing and transport of a pollutant is dominated byone-dimensional constant width flow, with velocities and mixing occurring primarily in thehorizontal plane only (simulating shallow water conditions) and only along the axis of the inlet. A1.5-D channel is essentially an extension to the 1-D case such that the channel width may vary butit is assumed that the flow in the cross-channel direction, y, is negligible. The 2-D case assumesboth flow and transport in both the x and y directions to be significant.Chapter 1 - Introduction5Within the context of this thesis, research has been conducted on one-dimensional situations.More specifically and particularly applicable to the British Columbia coast, this corresponds to thestudy of transport in long narrow inlets.Domain of Study CaseOf the influencing factors listed in Section 1.1, a few have been selected and studied in thisthesis. For the one-dimensional case the instantaneous release of pollutant for various channelgeometries and at various stages of the tidal cycle have been studied. In particular, thecharacteristics studied within this thesis include,1) 1-D channel2) instantaneous point-source release of pollutant3) a neutrally buoyant pollutant4) no return flow of pollutant5) no hydraulic forcing from additional inlets to the channelTo simplify the situation, return flow of pollutant on the tides was neglected. Also, the effects ofriver inflow are not studied. These do not, however, detract from the usefulness of the results.For example, the results would be directly applicable to areas where longshore currents account forlittle or no return flow and where runoff is minor. The instantaneous release would mimic apollutant spill or a periodic short-duration release from an industrial outflow pipe. The assumptionof neutral buoyancy is valid for pollutants that readily mix into the water column because of similardensity or those that are hazardous at very minor concentrations (such as PCB's) and are perhapssuspended in the water column.ObjectivesThe objective of the thesis is to improve the ability to predict the fate of a pollutant released in anarrow channel under the conditions described above. This has been achieved by developingsuitable numerical models which predict pollutant dispersion scenarios for various channelconditions. Numerical models have been developed which incorporate hydrodynamic componentsChapter 1 - Introduction6(describing the movement of water) and dispersion components (describing the movement of thepollutant in the absence of a current) for one-dimensional situations.Section 1.3 outlines previous work related to tidally forced pollutant transport and, to someextent, flushing processes in channels and harbours. Chapter 2 presents the theory behind thetransport mechanisms at work, and develops closed-form solutions for the hydrodynamiccomponents of the numerical model. Chapters 3 and 4 respectively discuss the numerical andphysical models undertaken. Chapter 5 compares the results obtained from Chapters 3 and 4 anddiscusses the influence of various parameters important to transport. Chapter 6 presents anoverview of the research conducted and provides some considerations for further study related tothis topic.1.3 Review of Pollutant TransportTidal Motions In the Pacific Northwest, tidal forcing is a primary mechanism for the disposal of effluent fromcoastal channels. A typical high tide exchanges a significant portion of the volume containedwithin the channel at slack tide (tidal volume). This is also related to the tidal prism - theproportion of the new water entering a channel compared to its tidal volume.In this region, tides are primarily diurnal mainly mixed (i.e. unequal tidal highs and lows with a12 hour period). From a modelling viewpoint, it is difficult to define a characteristic tidal cycle.To determine a characteristic tidal flushing sequence one might be obliged to model a full lunar tidalcycle. This of course would represent a period much greater than the residence time of most smallchannels. In Nece's model studies of tidal transport in harbours (1984), three tidal ranges wereused to characterize transport variations:1) neap tides2) mean tides3) spring tidesChapter 1 - Introduction7Another consideration relates to seasonal variations in local tributaries, i.e. large spring and autumnfluctuations in river outflow, which will influence the residence time of effluents.For a two-dimensional case, Nece (1984) describes the importance of gyres on the flood tide andhow they will influence the mixing process. Here, it seems necessary to assess the efficiency ofturbulent mixing processes because the effectiveness of tidal flushing is dependent on the rate atwhich the tidal flow mixes with the channel water. This, however, is not an easy task and someresearchers prefer to use a more general approach.Van de Kreeke (1983) explores harbour-related tidal exchange through the assumption that"currents induced by the wind, tide and density ... cause the waters in the basin to be completelymixed at all times". Within this context he describes residence time and flushing time for discreteand continuous dye injections in tidal prism models. He concluded that the amount of timerequired to expel the majority of the pollutant (residence time) for this completely mixed case is afunction of two dimensionless parameters:1) the fraction of the flood volume that is new water (i.e. tidal prism)2) the phase of the tide at which discrete injections are releasedThis approach (i.e. assuming completely mixed conditions) neglects the temporal importance ofmixing to the exchange process, but the tidal prism method does appear to provide a reasonablefirst approximation for estimating residence times.Wind and Wave EffectsThe effects of wind on mixing within a semi-enclosed body of water are much more limited thanin open water conditions. Within the context of this thesis, it is assumed that wind-driven wavesare of no consequence due to the sheltered nature of the case described. The wind forcing will ofcourse control the motions of the surface layer and of floatables such as oil slicks. Schwartz andImberger (1984) concluded that wind strength (as opposed to basin geometry relative to windChapter 1 - Introduction8direction) is the dominant factor controlling surface mixing. Wind shear will dictate the degree ofstratification in channels with fresh water sources. Furthermore, wind shear may cause tilting ofstratified interfaces and seiching effects when fluctuations in the wind shear forces occur.By nature, waves travelling through the seaward boundary do not frequently affect mixingfurther upchannel. Effectively, the source of wave action to be considered is that associated withvessel motions within the channel. These forcings will not affect flushing but they may inducespatial mixing. Depending on the channel depth relative to wave length, a wave's orbital motionmay cause resuspension of benthic material or simply aid in suspending non-buoyant material untilit is flushed out.Environmental Effects Related to Pollutant Transport in HarboursIn Washington State, regulations require that in the design of a channel, consideration is given tothe effects on juvenile salmon that migrate along the shorelines' shallow reaches (Nece, 1972).Within the confines of some local harbours, pollutant concentration levels reached toxic levels forthe vulnerable salmon fry. Other single entrance designs occupied the entire breadth of shallowzones and forced salmon fry to enter deep water areas where they fall prey to larger fish. In orderto prevent this, multi-entrance channels and detached breakwaters are now favoured over singleentrance enclosures to increase overall pollutant transport rates.In Australia, regulations require that potential impact to the water quality from such dredgedharbours and channels be assessed prior to construction (Schwartz and Imberger, 1986). A localconcern is the accumulation of seagrasses in these enclosures. Large amounts of organic materialmay grow in the semi-stagnant waters within channels, drastically reducing oxygen levels (i.e.imposing a high B.O.D.) as the material decomposes. Fish kills are related to the presence of these`dead zones' where mixing is poor and pollutant concentrations and low oxygen levels reach lethallevels. These are most commonly found at the end of long reaches. It is therefore important thatchannel and harbour designs avoid creating such zones.Chapter 1 - Introduction9It should be noted that regular tidal exchange within harbours and channels is beneficial to thesurrounding waters by allowing local organisms to adjust to the relatively dilute effluent. Ifconcentrated 'slugs' of polluted water exit the channel or harbour on an intermittent basis then themortality rates of local species would be expected to increase.Physical Models and Design The question remains as to whether it is more effective to model pollutant mixing and transportby means of numerical models or physical models. In reviewing the literature, such as Nece(1984), it seems that economical computer models that suitably incorporate velocity fields withdiffusive effects, did not exist at that time. Numerical models generally require well-definedboundary conditions in order to determine the importance of various input parameters. Physicalmodels are well suited to determine the relative importance of these parameters since each variablecan be altered to study its influence on the system.Nece (1984) describes further physical models of harbours, with which he studied the influenceof a number of aspects such as:1) entrance orientations, widths, and locations,2) various channel geometries (breadth:length ratios)3) the effect of rounding basin corners.Nece's primary interest in conducting these tests was to determine those variables important inmaximizing the overall tidal exchange to the large tidal range in the Pacific Northwest. Modelstudies were then compared to field studies for existing small-boat harbours, and results werepresented in terms of spatially averaged tidal exchange coefficients (i.e. the degree of flushingrelative to harbour width and length).Field StudiesSchwartz and Imberger (1986) describe, in detail, field studies conducted within a large marina(moorage for 1000 vessels). Using an extensive field program, they quantified a vast array ofChapter 1 - Introduction10physical parameters. Flushing rates of Rhodamine dye were compared to those described by thetidal prism method to determine the relative importance of the tide (with tidal range ofapproximately 1 m). Stratification within the layers was a prominent feature of this study and wasaffected to a large degree by wind shear conditions. During some ebb tides, the presence ofinternal supercritical flow within the constriction of the channel entrance was noted.In this particular case, it was concluded that the dominant flushing mechanism was infiltration offresh groundwater that aided in the mixing of the deeper stratified layers, thereby increasing theeffective flushing rates. Bienfang (1980) noted a similar occurrence where groundwater infiltrationincreased flushing rates by a factor of six to ten. These cases show that seemingly obscure factorsmay turn out to be dominant for specific sites.The field studies of Nece and Richey (1972) dealt with tracking miniature drogues that floated atdifferent subsurface heights. Other larger surface drogues were observed as they drifted aroundthe channel during ebb tides. These results (shown as spatial pathlines) were compared to smallscale physical models.The cited works indicate that there are various methods of approaching the subject of tidallyforced pollutant transport. There is no consensus as to the best means of predicting residencetimes of conservative pollutants (i.e. those that do not decay) or other associated parameters.While some authors choose to assume that mixing occurs instantaneously, others dedicate entirefield studies to the subject.Chapter 1 - IntroductionChapter 2 - TheoryAs described earlier, the problem of numerically modelling pollutant transport in tidal situationscan be divided into two components:1) Hydrodynamic2) Advection/Diffusion (or Dispersion) TransportThe governing hydrodynamic equations shall be developed separately from the pollutant transportequations and give rise to the flow velocities and water elevations within the channel. However,the pollutant transport equations are governed by fluid motion and hence, require thehydrodynamic solutions within their own development.Figure 2.1 provides definition sketches describing the notation used herein. A channel ofconstant average depth, d, and length, X, is closed at one end, x = X. The channel has a width, w,and pollutant enters the channel at position Xp. The tidal wave is described using an amplitude A,length L, celerity c, and wave period T (= L/c). The surface elevation is denoted asThe assumption of one-dimensional flow considerably simplifies the pollutant transport problem.As mentioned previously, variable-width one-dimensional channels provide a way of predictingflows by assuming that flows in one direction (commonly the x direction) are of orders ofmagnitude larger than those in other directions (the y and z directions). Much research has beenconducted on fluid flow and transport mechanisms in two- and three-dimensional situations andthese involve complex solutions that are still based on simplified boundary conditions. As theboundary conditions become more structured and case specific, the equations to the hydrodynamicand transport components become increasingly difficult to solve.11Chapter 2 - Theory122.1 Governing EquationsFigure 2.1 provides definition sketches of the situation being considered.Hydrodynamic EquationsIncompressible fluid flow in three dimensions is described by the Navier Stokes equations andthe continuity equation given below (e.g. Sarpkaya and Isaacson, 1981).Equations of motionauau au au aU —a. + V — W — — (p + pgz) + DV2U [2.1]at ay az p axavav av av a + U — + V — + W — = — — (p + pgz) + DV 2V [2.2]at ax ay az p ayaw aw aw aw 1 a— + U + V + W — (p + pgz) + oV2W [2.3]at ax ay az p aza2u a2u a2uwhere V2 U= ax2 ay2 az2Continuity Equationau av aw+ — + = uax ay az[2.4]where U, V, and W are the velocity components in the x, y, and z directions, respectively; p is thefluid density; g is the gravitational constant; p is the fluid pressure; and v is the kinematic viscosityof the fluid.Assumptions A number of assumptions may be made in order to reduce the governing equations to a simplerform so as to apply to a one-dimensional channel. First, the tidal wave amplitude, A, is assumedto be sufficiently small compared to depth so that nonlinear terms in the governing equations can beChapter 2 - Theory[2.8]13neglected. Secondly, the tidal wave length, L, is assumed to be significantly greater than the waterdepth, d, so that the vertical particle accelerations may be neglected. In addition, energydissipation terms, including viscous effects giving rise to bottom frictional resistance and windstress effects, are neglected. Finally, the channel is assumed to have a constant average depth anda constant average cross-sectional area, Ar.Thus the assumptions made are that:1) A « L, d2) L>>d3) no friction4) constant Ar, i.e. V = 0By applying these assumptions to the equations of motion and omitting the nonlinear terms, [2.1]to [2.4] reduce to:au 1 ap- —at = p axaw = - —1 -a-- (p + pgz)at p azau aw 0ax + az - -These represent the x and z momentum equations and the continuity equation, respectively.DispersionThe conservation of mass for the pollutant gives rise to the following governing equation:d (ArC) a (ArC) a (UArC) a (._ + .., - ax DAr ac + KdC +dt at ax ax qwhere C is the depth-averaged and width-averaged pollutant concentration, A r is the cross-sectional area of the channel, D is an eddy diffusion coefficient, Kd is a decay coefficientapplicable to non-conservative pollutants, and q are source/sink terms (Koutitas, 1970).Chapter 2 - Theoryat z = 0at z = 0W = 0 at z = -danw . —atp = 014AssumptionsIn the present context, the sink/source terms are neglected so that the pollutant is introducedinstantly at time t = 0. Furthermore, the diffusion coefficient, D is assumed independent oflocation along the channel, and finally it is assumed that the pollutant does not grow or decay (i.e.it is conservative).By neglecting these terms, [2.8] reduces to:a(Arc) + a(UArC) _ D a (A acr )at ax ax - ax [2.9]This allows for variations of cross-sectional area, velocity and concentration with both x and t.2.2 Initial and Boundary ConditionsBefore solutions may be found for the hydrodynamic and pollutant transport components, theboundary conditions must be stated. In describing these, reference is made to Fig. 2.1.HydrodynamicNo initial conditions are necessary for the hydrodynamic components since the fluctuations ofTi(x,t) and U(x,t) are considered to be in a periodic steady state. The relevant boundary conditionsin the vertical direction are:These are respectively the kinematic and dynamic free surface boundary conditions at the stillwater level, and the kinematic boundary condition at the seabed. The free surface boundaryconditions are applied at z = 0 rather than z = 11 on the basis of the linearization entailed in the longwave theory deriving from the assumption A << d.Chapter 2 - Theory15Further boundary conditions required at the seaward and landward extents of the channel (x = 0and X, respectively). Open-water tidal fluctuations are assumed sinusoidal and may be describedby:is = A cos (wt) for t 0 [2.13]where TN is the water elevation outside the channel, w is the tidal angular frequency and A is thetidal amplitude. The water level at the seaward boundary is taken to be equal to the open-watertidal water level:TI (0,t) = 11 s [2.14]At the landward boundary (x = X) the fluid velocity is taken to be zero:U(X,t) = 0 [2.15]Dispersion The initial condition for an instantaneous point discharge of pollutant at location x p in a channelwith a zero ambient concentration is:0 for x # xpC(x,0) = [2.16]{ Co 8(x) for x = xpHerein 8(x) is the Dirac delta function and Co is given as:M Co =PAr[2.17]where M is the mass of the pollutant discharged. The Dirac delta function is defined such that:+00f 6(x) dx = 1 [2.18]Chapter 2 - Theory16At the initial time, the concentration corresponds to an impulse of zero width and total pollutantmass, M. The boundary condition at the landward boundary is such that the pollutantconcentration gradient here is zero, so that the pollutant cannot be transported through theboundary walls. Thus:ac— = oanat x = 2t, [2.19]where n is the distance normal to the solid boundary (i.e. along x).At the seaward boundary (x = 0), the boundary condition is such that the pollutant flux gradienthere is zero (and inflow concentrations are zero) (Koutitas, 1970).acD-- - UC < 0axat x = 0 [2.20]Furthermore, the concentration here is equal to the open water ambient concentration, so that thereis no diffusion of pollutant through the seaward boundary toward the channel. As alreadyindicated, this is taken to be zero, so that:C(x,t) = 0 for x<_0 [2.21]2.3 Development of Hydrodynamic ComponentThe hydrodynamic component developed herein provides the current velocity, U(x,t) and waterlevel, ri(x,t) at various positions and times within the channel in terms of the following parameters:1) average water depth, d2) channel length, k3) tidal amplitude, A4) tidal period, TA closed-form solution to the hydrodynamic equations for the case of interest may be developedusing linearized long wave theory (e.g. Sarpkaya and Isaacson, 1981). On the basis of theChapter 2 - Theory17previous assumptions, the horizontal particle acceleration is invariant with depth and by virtue of[2.6] the fluid pressure is hydrostatic:p = pg (i-z)Substituting this into [2.5] provides:au ailat + g Tx =0[2.22][2.23]This implies that horizontal velocity is uniform over depth, U = U(x,t). By integrating [2.7] fromz = -d to z = T1, applying the boundary conditions [2.10] to [2.12] and assuming that it << d, weobtain:all + d au = oat ax [2.24][2.23] and [2.24] can be combined to form the following second order equations given in terms ofa celerity, c.a211 2 a211 = 0- c2 [2.25]at2 ax2where:a2u a2u - C2 d — oat2 ax2 [2.26] c = I g d [2.27]Equation [2.25] is commonly referred to as the wave equation. The general solutions to [2.25] and[2.26] are of the form:T1 = fl (x - ct) + f2 (x + ct) [2.28]U = f3 (x - ct) + fg. (x + ct) [2.29]Chapter 2 - Theory18which correspond to the superposition of two steady wave forms propagating at a speed c, one inthe positive x-direction, and one in the negative x-direction. If both 11(x,0) or 'n(0,t) are known,then the subsequent motions for all times can be determined. Because the boundary condition,[2.14] is harmonic, [2.28] and [2.29] can be expanded to:11 = al cos (kx - cot) + a2 cos (kx + cot) + a3 sin (kx - cot) + a4 sin (kx + cot) [2.30]U = c1 cos (kx - cot) + c2 cos (kx + cot) + c3 sin (kx - cot) + c4 sin (kx + cot) [2.31]where k = 27c/L, co = 27t/T, and an and e n are constants which may be determined for the landwardand seaward boundary conditions.The derivative of [2.30] with respect to x and the derivative of [2.31] with respect to t can nowbe substituted into [2.23]. This produces a relationship between an and cn.c 1 = al a, c2 = -a2 a, c3= a3 a, c4 = - a4 a [2.32]where a = gk/co = g/c. By applying the landward and seaward boundary conditions [2.14] and[2.15], expansions for an may be developed as:al = a2 = A [2.33]A ,, , ,a3 = a4 = -2-, tan .1(A..) [2.34]By substituting [2.33] and [2.34] into [2.30] and [2.31] respectively, the closed-form solutionmay finally be developed as:1 (x,t) = A cos (cot) [ cos (kx) + tan (kA) sin (kx) ] [2.35]U (x,t) = -g- sin (cot) [ sin (kx) - tan (k?) cos (kx) ] [2.36] Chapter 2 - Theory19The closed-form solutions for ri (x,t) and U (x,t) obtained here have a considerable advantagein not contributing to any approximation errors which arise in the finite difference methodssubsequently employed in solving for the advection/diffusion components (described inChapter 3).2.4 Dimensional AnalysisFinally, it is useful to carry out a dimensional analysis for the concentration, C, in order toestablish those parameters which influence the problem at hand. The concentration may be taken todepend on the following independent variables which define the problem:C (x,t) = f ( x, xp , t, to, d, X, T, A, g, D, Co ) [2.37]A dimensional analysis then provides the following relationship between dimensionlessparameters:C(x,t) f ( x t xp to d A T2g D2 \Co 2 ' d' A'gd3 ' [2.38]where C/Co is the normalized concentration; x/X and UT are dimensionless measures of theposition and time respectively; x p/X and to/T are the position and instant of the initial pollutantdischarge; dA is the ratio of the water depth to channel's length; A/d is the ratio of tidal amplitudeto water depth; T2g is a measure of the tidal wave length to channel length ratio; and D 2/gd3 is adimensionless measure of the diffusion coefficient.Chapter 2 - TheoryChapter 3 - Description of Numerical ModellingAn initial step in the research was to develop a one-dimensional numerical model to assist in abasic understanding of the physical processes at work. The model is required to predict pollutanttransport within the tidally-forced currents of a one-dimensional inlet of uniform depth and width,with no river forcing or stratification. This case is depicted in Fig. 2.1.The numerical model is based upon the closed-form solution for the hydrodynamic componentand finite difference approximations for the dispersion component. Previous one-dimensionalmodels were not used as a basis for the research since it was intended that the process ofdeveloping such models would aid in understanding the fundamentals of tidally-forced pollutanttransport.The numerical model was written in FORTRAN and run on a VAXstation 3200 computer.Numerical results from these models were visualized using the graphics capabilities of Macintoshand Sun sparcstations computers.Testing Procedure The numerical model was developed in stages in order to allow for an understanding of eachmechanism involved and in order to ensure the reliable development of the complete model.Firstly, the hydrodynamic component of the model was developed and tested. Secondly, withinthe pollutant transport component of the model, the advection and diffusion components weredeveloped separately in order to study the nuances of numerically modelling each of thesemechanisms. They exhibit different characteristics depending on the application of similar finitedifference methods.20Chapter 3 - Description of Numerical ModellingOnce the features of the component models were adequately understood and tested, a complete 21transport model was compiled and tested for a situation involving uniform velocity (i.e. no tides);and only then was it combined with the hydrodynamic component so as to model the entire tidally-forced situation. Thus the testing procedure that was adopted is as follows:1) hydrodynamic component [ numerical model # 1 ]2) i) advection component [ numerical model # 2 ]ii) diffusion component [ numerical model # 3 ]iii) advection & diffusion, i) + ii) [ numerical model # 4 ]3) coupled components, 1) & 2) [ numerical model # 5 ]The details of each test scheme are described in Table 3.1 to 3.3. Further explanation of thevariables are given in the following sections.3.1 Hydrodynamic ComponentThis section of the model was based on the closed-form solution, [2.36] to [2.37], developed inSection 2.3. Since the closed-form solution is exact for the given boundary conditions, nosensitivity analyses were needed such as those conducted for the following advection and diffusioncomponents which use finite difference techniques. The output from the hydrodynamic modelprovides the velocity U(x,t) and the surface elevation n(x,t) at a series of positions and instancesand was displayed in two dimensional arrays.3.2 Advection/Diffusion ComponentThe governing dispersion equation [2.9] for one-dimensional variable-area pollutant transportmay be simplified for cases where the width is uniform and elevation changes are small relative tothe depth. Equation [2.9] then simplifies to:ac + a (UC) _ 2cat ax ax2 [3.1]Chapter 3 - Description of Numerical ModellingFurthermore, if the velocity, U, is taken as uniform along the channel, then [3.1] gives:ac ac a2c— + u — = Dat ax ax2 [3.2]Equation [3.2] differs from [2.9] in that the hydrodynamic component is omitted, resulting inconstant volume and velocity conditions within the test section.Finite Difference MethodThe finite difference method (FDM) is one of the primary means of transforming the partialdifferential equations (PDE's) governing fluid flow to a predictable discrete system of linearequations. These equations may then be solved numerically at a series of discrete time steps. Theaccuracy of this method depends greatly on the size of the finite differences implemented. As thesizes of the discrete intervals grow, so does the discrete (or approximation) error. Conversely, asspatial and temporal spacings approach zero the finite difference scheme approaches the continuoussystem of the PDE's.The schemes used herein are all explicit finite difference schemes, such that any value at anadvanced time t + At may be calculated explicitly from the known solution at time t. Explicitschemes generally suffer from the phenomenon of truncation error, suggesting that the pollutant isartificially advected or dispersed to a greater extent than is implied by the governing equations.The other type of scheme commonly employed for reducing PDE's is known as an implicitscheme. Implicit schemes entail the solution of unknowns within large algebraic matrices, andgenerally require considerable computing power.As stated earlier, exact solutions to the governing equations for concentration are not possible,so that finite difference techniques are used to approximate the required results. The accuracy ofsuch results depends upon two factors: the numerical scheme employed and the relative values ofnumerical parameters such as the time step size.22Chapter 3 - Description of Numerical ModellingConsider a function f(x,t) at discrete values of x and t with intervals of Ax and At, respectively.In the simplest schemes based on forward, backward, and central differencing, the first-orderderivative df/dx is approximated as follows:Forward Differencingdf f(x+Ax,t) - f(x,t) dx — Ax [3.3]Backward Differencingdf f(x,t) - f(x-Ax,t) dx — AxCentral Differencingdf f(x+Ax,t) - f(x-Ax,t) dx — 2AxA second-order derivative, d2f/dx2 , can be approximated in a similar manner. In particular, thecorresponding central difference is:d2f f(x+Ax,t) - 2f(x,t) + f(x-Ax,t) By expanding f(x+Ax,t) as a Taylor series in terms of f(x), df/dx, etc., it can be shown that theforward and backward difference schemes listed as [3.3] and [3.4] are first-order accurate, withapproximation errors denoted as O(Ax). The central difference schemes [3.5] and [3.6] aresecond-order accurate, with approximation errors denoted as O(Ax 2). For [3.3] and [3.4], thisimplies that as the interval Ax is halved the associated errors are halved, whereas those errorsassociated with [3.5] and [3.6], which are second-order accurate, are quartered.The approximation (or truncation) errors for the various pollutant transport models may bedetermined through Taylor series expansions of the FDM operators contained within the transport23[3.4][3.5]dx2 — Ax2 [3.6]Chapter 3 - Description of Numerical Modellingequations. Different combinations of finite difference operators create different truncation errorsof varying magnitudes. First-order truncation errors affect the phase of the pollutant peak as itadvects along the channel's length and may be referred to as artificial advection. Second-ordertruncation errors affect the diffusion of the pollutant and are referred to as artificial diffusion.Higher order truncation errors are of less importance than first- and second-order errors and aremore difficult to associate to physical phenomena such as advection or diffusion. It should benoted that zeroth order truncations are possible and affect the bias of the results. If this occurs forthe cases in question, this would suggest that the mass of pollutant is not being conserved.Figure 3.1 shows the finite difference scheme employed for the combined model. It is possibleto produce higher accuracy terms for forward and backward differences, but as the accuracyincreases so does the need for information from cells that are further away from the grid point inquestion.Initial ConditionsIn testing the various finite difference models developed and in the associated laboratory tests,the intent is to apply an initial condition corresponding to the instantaneous discharge from a pointsource. However, in the laboratory, a finite amount of time (-- 30 sec) is required to discharge thepollutant during which spreading occurred and, subsequently, there is a finite width to the pollutantcloud (-- 40 - 50cm) at t = 0. Moreover, numerical methods such as FDM have trouble simulatingextreme gradients such as one-point spikes and Dirac functions, and rely on concentrationinformation to be presented less abruptly over several grid points thus producing shallowerconcentration gradients.For this reason, several different initial pollutant shapes were tested to determine the smallestnumber of points needed to accurately define the input profile. Results have been obtained for four24Chapter 3 - Description of Numerical Modellingdifferent initial concentration profiles (i.e. at t = 0) as sketched in Fig. 3.2. These correspond to:1) one-point triangle2) three-point triangle3) five-point triangle4) three-point rectangleOf these shapes, the one-point triangle and the three-point triangle were rejected owing to thedifficulty in modelling these numerically using first-order central differencing. This is due to theirsevere initial slopes in the concentration data. Even though the purpose of the research is to modelinstantaneous effluent releases from point sources (which would be represented by a very narrowrectangular block of pollutant within the channel), the five-point triangular solution shapeadequately resembled the effluent input in the physical model and still remained a manageable sizein the initial tests. Within the advection, diffusion, and advection/diffusion models a five-pointtriangular concentration profile is therefore used.For the final numerical model, however, a larger point triangle was eventually used based on agrid size that was small enough model local variations concentration and also the initial pollutantcloud width within the physical model. Based on a grid size of 5 cm and estimating the averagepollutant cloud width at approximately 45 cm, a nine-point triangle was finally chosen. The largernumber of points reduces the possibility of problems within the numerical scheme.3.2.1 AdvectionAdvection is defined as "...transport [of a pollutant] by an imposed current system..." (Fischer etal, 1982). In order to examine the case of advection alone, diffusion is considered not to occur.The corresponding governing equation is essentially [3.2] with the diffusion coefficient set to zerofor the case of uniform velocity along the channel, and this is simplified to:acac— + u — = 0at ax[3.7]25Chapter 3 - Description of Numerical ModellingA general solution to [3.7] is:C = f (x - Ut)which corresponds to the pollutant concentration maintaining a constant spatial profile whichmoves at a speed U.Finite Difference MethodSince the terms used in [3.7] are first derivatives, the advective process can conveniently bemodelled using a forward-in-time, backward-in-space (FTBS) finite difference operator.Equations [3.3] and [3.4] are applied to the first derivatives in [3.7]. This provides:C (x,t+At) - C (x,t)C (x t) - C (x-Ax,t) — - U 'At AxThe concentration at the advanced time t +At is therefore expressed as:C (x,t+At) = C (x,t) - UAx [C (x,t) - C (x-Ax,t) ]AtUAt is known as the Courant number, Cr. Equation [3.9] may be rearranged as:AxC (x,t+At) = ( 1 - Cr) C (x,t) + C r C (x-Ax,t) [3.10]Boundary ConditionsIn carrying out tests on advection with a uniform velocity, the boundary conditions with respectto the hydrodynamic component are absent. A simple periodic boundary condition was thenapplied to the advection equation such that any pollutant which moved across the right-handboundary of the domain in any time step, reappeared at the left-hand boundary. This implies thatafter one cycle, of period X/U, the pollutant should be in its initial position with an unchangedconcentration profile.26[3.8][3.9]Chapter 3 - Description of Numerical ModellingTruncation ErrorFrom the Taylor series expansion of [3.7] the associated truncation error, E is:a cE = 1-2 Ax U (1 - Cr) ax2 27a 3 C + I Ax2 U (1 - Cr)2 + (higher order terms)6 ax3 [3.11] This suggests that artificial diffusion (from the a 2C/ax2 term) and higher order errors exists forthose results where Cr 1. On the other hand, for Cr = 1, the FDM scheme is exact and all higherorder terms are eliminated.Numerical Tests - Effects of CrSince it is known that the Courant number influences the solution for the concentration (e.g.Abbott and Basco, 1989), the effects of Cr were examined by carrying out tests using a five-pointtriangular dye slug to represent the initial concentration and width with a channel length, X = 600m, and velocity, U = 0.75 m/s and various values of Cr. Table 3.1 shows the testing programadopted for this case.By varying the Courant number, three different categories of results are produced.1) Cr « 1.0: solutions are heavily damped and insufficiently accurate.2) Cr 1.0: solutions are identical or reasonably close to the exact solution.3) Cr > 1.0: solutions are unacceptable or completely unstable due to strongoscillations.Corresponding concentration profiles at the initial times and one cycle later when t = X/U arecompared in Fig. 3.3 for 4 values of Courant number, Cr = 0.5, 0.75, 1.0, 1.5. As seen inFigures 3.3 (a) and (b) (Cr = 0.50 and 0.75, respectively), the peak amplitudes are lowered andthe shape of the dye slug has changed. Note that no phase lag has occurred in either case but theheavy damping makes the results too inaccurate to be of much use. Figure 3.3 (c) (C r = 1.0)Chapter 3 - Description of Numerical Modellingshows that the solution surface of the propagating dye slug is identical to the initial shape. Thiscase relates to [3.10] reducing to the following:C (x,t+At) = C (x-Ax,t) [3.12]which may be seen to be equivalent to:C = f (x - Ut)when Cr = 1.For Cr = 1.5, shown in Figure 3.3 (d), the returning solution surface bears little resemblance tothe initial shape due to strong oscillations. The instability of this case increases as time progressesand is clearly unstable.From these tests, it is clear that for a case of pure advection the Courant number should be heldas close to 1.0 as possible without exceeding this limit If C r is likely to vary for a given situationwithin a predetermined range, it is better for C r to err on the low side and have numerically dampedresults than to have results that are unstable.Effect of grid size and time stepCertain aspects of the selection of the variables within the Courant number (i.e. U, Ax, and At)must be mentioned. Assuming that the velocity, U is specified a priori by hydrodynamicconsiderations, one must select the ratios of Ax and At to provide suitable values of Cr. The valuesof Ax and At should be selected to obtain results that are sufficiently accurate, but at the same timeare efficient in terms of computing time and storage. It will be seen in subsequent sections thatwhen diffusion is considered, the ability to vary Ax and At is greatly reduced. In a more specificcase, tidal flow velocities oscillate between zero and a maximum value described in Section 2.3.This implies that the Courant number will also vary accordingly, assuming Ax and At are constant.Figure 3.4 shows the test case (X, = 600 m, U = 0.75 m/s) where the Courant number is heldconstant at 0.75 while the grid sizes (Ax and At) are varied. Although C r is held at a less than ideal28Chapter 3 - Description of Numerical Modellinglevel, the FDM used returns solution shapes which approach the exact shape as the grid sizes arereduced. This implies that for values of C r less than one, [3.8] is convergent. Stated otherwise,as spatial and temporal spacings approach zero, the FDM terms approach the corresponding termsof the PDE:29AAtC AC+ U Ax =>aC ac+ u at ax as Ax, At => 0 [3.13]Interestingly enough, if the same procedure is applied to a case with a Courant number greaterthan one, the same reduction in grid size results in a decrease in stability for [3.8]. This caseindicates a divergent nature.3.2.2 DiffusionThe case of diffusion alone, without advection, may be examined by setting the velocity, U equalto zero in [3.2], giving rise to the simplest governing equation for diffusion:ac= D a2cat axe [3.14]Exact SolutionAn exact solution exists for this situation against which the results of the FDM may be compared.For this case, the pollutant released at t = 0 has an infinitely small width and an infinitely highconcentration. The area under this initial concentration profile is described by the Dirac function.For the case of no boundaries, C (±00, t) = 0 and for t > 0, the closed-form solution (e.g. Fischer,1981) is:C (x,t) = exp (- (x4-DxPt )2 )(Ai 41‘47tDt[3.15]where M is the mass of the discharged pollutant at t = 0, and xp is the position of initial dischargeof the pollutant. Although the initial concentration profile used in the numerical models (with aChapter 3 - Description of Numerical Modellingfinite width) differ somewhat from the ideal initial condition used above, it is anticipated that this 30should not affect the exact solution which results after a sufficient duration.Finite Difference MethodIn contrast to the advective component which was modelled using a forward-in-time, forward-in-space operator, the first and second derivatives within the diffusion equation are most easilymodelled using a forward-in-time, centred-in-space (FTCS) finite difference operator. Theboundary elements are handled slightly differently and are discussed below. The FTCS operatorimplies that spatial information on both sides of a given grid point is needed in order to determinethe concentration at subsequent time steps. Equation [3.14] is therefore approximated by using[3.3] and [3.6] and provides:C (x,t+At) - C (x,t) C (x+Ax,t) - 2C (x,t) + C (x-Ax,t) — D [3.16]At Ax2Again, solving for the concentration at the advanced time C (x,t+At):DAtC (x,t+At) — 2 [ C (x+Ax,t) - 2C (x,t) + C (x-Ax,t) ] + C (x,t)60( [3.17]By denoting DAt/Ax2 as the dimensionless parameter y, [3.17] may be written in the form:C (x,t+At) = y C (x+Ax,t) + (1 - 2y) C (x,t) + y C (x-Ax,t) [3.18]Equation [3.18] indicates that the concentration at a given location at the advanced time t+At,depends on the concentrations at that location and the two adjacent grid points, all at the previoustime step, t.Boundary ConditionsA numerical model using a five point triangular pollutant slug is again employed to test this case.One set of boundary conditions was used to model this channel section, i.e. closed right and left-Chapter 3 - Description of Numerical Modellinghand boundaries:aC/an = 0 at x = 0 and x = X [2.19]A closed boundary condition is achieved numerically by placing an additional numerical cellbeyond the model's end, and equating the concentration in this cell to that in the adjacent (true) endcell:C (X+Ax,t) = C (A,t) [3.19]This technique, referred to as mirror imaging, allows the FDM to be applied to the end cellwithout modification and establishes a zero concentration gradient at the end cell, thereby satisfying[2.19]. To test the integrity of the closed boundaries the amount of pollutant injected into to thechannel should remain constant. The mirror imaging technique provides this condition andsatisfies the conservation of mass for the pollutant.Truncation ErrorFrom the Taylor series expansion of [3.14], the associated truncation error, E, is:1 a4C E = D Ax2 ( 6— - y) + (higher order terms)ax4 [3.20]Since the first term is of fourth order, this suggests that very little error is associated with [3.16]and that no artificial advection or diffusion (being first and second-order errors, respectively) existwithin the scheme. Nevertheless, [3.20] may still be optimized for least error corresponding toy= 1/6. For this case, [3.14] becomes fourth order accurate in space and second-order accurate intime, i.e. O(Ax4 , At2).Numerical Tests - Effects of yThe parameter y has similar tendencies as the Courant number (Section 3.2.1) in causing theFDM to become unstable beyond a threshold value - in this case 7 = 1/2 (Abbott and Basco, 1989).This value of y, however does not provide the most accurate results for the diffusion case as31Chapter 3 - Description of Numerical Modellingshown above. The pure diffusion case is therefore bounded by one y value and most accurate foranother. The range of y values between 0 and 1/2 is still third order accurate in any case and assuch is extremely precise when compared to the other schemes tested.For the case of 7 = 1/6 [3.18] becomes:1C (x,t+At) = —6 C (x+Ax,t) + g4 C (x,t) + 1 C (x-Ax,t)6 [3.21]By varying y, three categories of results are produced:results are fourth order accurate.solutions are stable and are third order accurate.solutions are unconditionally unstable.Table 3.2 shows the testing procedure followed for the case of pure diffusion.Figures 3.5 display the results of diffusing the five-point slug over time using various 7 values.Although the results cannot be compared directly due to the use of various diffusion coefficients,they do show the stability of the results for a common time step and grid size. Figures 3.5 (a), (b),and (c) for 0 < 7 < 0.5 show an increasing amount of diffusion as y increases. Figure 3.5 (d), for7 = 0.55, displays increasingly unstable results with time.3.2.3 Advection/DiffusionNow that the simpler cases corresponding to advection alone and diffusion alone have beenaddressed, the case of combined advection and diffusion is considered. Again, the governingequation, for the case of a uniform velocity is:ac ac rcjat ax ax2 +._. D — [3.1]32Chapter 3 - Description of Numerical ModellingExact SolutionTo illustrate this case, a test series was run using FDM and compared to the exact solution(Fischer, 1981).C (x,t) = , M exp (- E u(x-xP ) 12 )4DtAl 4nDt[3.22]Finite Difference MethodThe FDM employed here uses a forward-time, centred-in-space (FTCS) operator in a slightlyaltered form compared to that used in the Sections 3.2.1 and 3.2.2. In this situation, wherediffusion plays a role in transporting the pollutant ahead of the advected slug, central differencingmay be employed to increase the accuracy of the advection term to second-order, O(Ax 2). This wasnot possible (or necessary) in a strictly advective case since the downstream cell contains noconcentration information. Equation [3.1] is therefore approximated by using [3.3], [3.5], and[3.6] and provides:C (x,t+At) - C (x,t) C— lj (x+Ax,t) - C (x-Ax,t) At 2Ax+ D C (x+Ax,t) - 2C (x,t) + C (x-Ax,t) Ax2 [3.23]Combining these terms and implementing the notation C r and y, the governing finite differenceequation is:C (x,t+At) = [7 — 0.5Cr] C (x+Ax,t) + [1 - 2y] C (x,t) + [y + 0.5Cr] C (x-Ax,t) [3.24]Stability range for y and CrThe stability range for y and C r is shown in Fig. 3.6 (after Leonard, 1979). Note that for y = 0the solution in unconditionally unstable. The curved upper boundary of the graph results from alinear stability analysis, LSA, of [3.24] and shows that Cr and y can each be within their individual33Chapter 3 - Description of Numerical Modellingthreshold values and still produce unstable results when they are both present within a FDMscheme.An additional dimensionless parameter shown in Fig. 3.6 is the Peclet number, Pe and may benow employed to describe the relative importance of advection and diffusion for a given set ofconditions (Abbott and Basco, 1989). This is defined as:P Cr U Axe = [3.25]If the upper limits of 7 and Cr are used (0.5 and 1.0, respectively) Pe attains a value of 2.0. Forvalues of Pe > 2.0 advection is of primary importance, whereas for P e < 2.0 diffusion is thedominant factor.Boundary ConditionsFor this test case of advection and diffusion with a uniform velocity, a five-point slug andperiodic boundary conditions are again employed. As in Section 3.2.1, any pollutant which movesacross the right-hand boundary of the domain in any time step reappears at the left-hand boundary.After one cycle, of period VU, the pollutant should return to its initial position with a concentrationprofile changed by the diffusion process.Truncation ErrorFrom the Taylor series expansion of [3.1] the associated truncation error, E, to the fourth orderterm is:a3cE At u2 at2 + 2 U Ax2 - 6 ) a,(3+ D Ax2 ( - y) a4c+ (higher order terms)6 ax4[3.26]34Chapter 3 - Description of Numerical ModellingThis suggests that the error associated with [3.26] contains artificial diffusion related to AtU 2arising from the advection term. To minimize truncation error, one should attempt to reduce At asmuch as possible, assuming that U is uncontrollable. It is also interesting to note that the third andfourth order component of E are eliminated for y = 1/6.Numerical Tests - Effects of y and Cr combinedFigure 3.7 shows the effect of varying Pe for the case of y= 1/6. Table 3.3 shows the testingprogram followed. Relating to the Leonard's diagram, Pe may vary between 0.0 andapproximately 3.7. One can see that the model is stable for Pe = 1.0 and 2.0 but for larger Pevalues (i.e. Pe = 4.0) the numerical results begin to oscillate. For the particular values of D, U, At,and Ax chosen, results most closely resemble the exact solution for Pe = 2.0. Figure 3.8 showsnumerical results for three y values (y= 1/12, 1/6, and 1/2) and various Pe values. Once more thestability of the results are seen to be governed by an upper limit relating to Fig. 3.6. It should bestressed that when dealing with advection and diffusion one is obliged to vary U, D, At, and Ax toobtain prescribed values of Cr and y.3.3 Combined Components and ExtensionsThe hydrodynamic, advective and diffusive components developed in sections 3.1 and 3.2 arenow combined, allowing both U and Ar to vary. In the corresponding model, which represents atidally forced one-dimensional channel, the velocity and channel cross-sectional area vary in timeas the tides rise and fall. For this situation no closed-form solution is known to exist.Equation DevelopmentThe governing equation for the corresponding model is:a aat (ArC) + ax (UArC) = D ax ax—a (Ar ac) [2.9]35Chapter 3 - Description of Numerical Modelling36For the case in question, the dispersion coefficient, D, is again assumed to be independent of itsposition along the channel, x/X.Equation [2.9] may be expanded as:ac aAr ac au aArAr + C + ArU + ArC + cuat at ax ax ax [3.27]D A a2c + D aAr acr — — —ax2 ax axExpressions for the derivatives of U and Ar appearing in [3.27] may be derived from the exactsolution [2.36] to [2.37]. Assuming that the channel is represented by width, w, the cross-sectional area is defined as: Ar(x,t) = [ d + I](x,t) w [3.28]Since d and w are constant, variations in the cross-sectional area, Ar, depend solely on the variations of the surface elevation, Expressions for the derivatives of U and A r are as follows:ax— Ar(x,t) = w — ri(x,t) = Akw cos (wt) [tan (kX) cos (kx) - sin (kx)]x a axa aatat Ar(x,t) = w — To,t) = - Mow sin (wt) [cos (kx) + tan (0) sin (kx)]ax gAk — U(x,t) = sin (cot) [cos (kx) + tan (kX) sin (kx)]c[3.29][3.30][3.31]Finite Difference MethodAgain, the FDM scheme used here is the FTCS for all grid points including the landward andseaward boundaries. Equation [3.3] is employed to model the forward-in-time component, ac/at,whereas [3.5] is used for the first-order central difference operator, aC/ax, and finally [3.6] isused to model the second-order central difference operator, a2C/ax2.Chapter 3 - Description of Numerical ModellingUsing the above abbreviations and implementing the FDM operators within [3.27], we obtain:C (x,t+At) - C (x,t) Ar(x,t) At + C (x,t) Art (x,t) + [3.33]C (x+Ax,t) - C (x-Ax,t)Ar(x,t) U(x,t) 2Ax + Ar(x,t) C (x,t) Ux(x,t)t) - 2C (x,t) + C (x-Ax,t)+ C (x,t) U(x,t) Arx (x,t) - D Ar(x,t) C (x+Ax,t) C (x+Ax,t) - C (x-Ax,t) - D Arx(x,t) 2Ax — 0where the subscripts x and t denote the corresponding derivatives. Making use of the parametersCr and y, [3.33] may be rewritten as:Art (x,t)C (x,t+At) = - Ar(x,t) C (x,t) AtrC - 2 [ C (x+Ax,t) - C (x-Ax,t) ]_ C (x,t) Ux(x,t) AtArx(x,t)- C (x,t) U(x,t) AtAr(x,t)[3.34]+ y [ C (x+Ax,t) - 2 C (x,t) + C (x-Ax,t) ]Arx (x,t) y Ax+ 2 Ar(x,t) [ C (x+Ax,t) - C (x-Ax,t) ]Gathering terms with respect to the cells x+Ax, x, and x-Ax, [3.34] may be expressed as:C (x,t+At) =A rx (x,t) y Ax)( y - 0.5Cr + 2 Ar(x,t) C (x+Ax,t) [3.35]( A rt(x,t) At Arx(x,t) U(x,t) At- Ar(x,t) + U x (x,t) At + Ar(x,t) + 2y ) C (x,t)A rx (x,t) y Ax2 Ar(x,t)+ (y + 0.5Cr - C (x-Ax,t)37Chapter 3 - Description of Numerical ModellingEquation [3.35] may be used to solve for the concentration at an advanced time C (x,t+At) usinga known solution at the present time t. Therefore, if the initial concentrations within the channelare known all subsequent concentration profiles may be determined.Boundary ConditionsAs described earlier, a nine-point slug is used to simulate the initial pollutant cloud from thephysical experiments. Two sets of boundary conditions were needed to model this channel sectionand are shown in Figure 3.1:1) Open left-hand boundary condition, C = 02) Closed right-hand boundary condition, DC/ax = 0The open left-hand boundary condition was achieved numerically by placing an additionalnumerical cell beyond the open boundary of the model and by setting the concentration here tozero. This provides a positive concentration gradient satisfying [2.22]. i.e.C (-Ax,t) = 0 [3.36]This would equate to a physical situation of a longshore current moving the pollutant away fromthe channel mouth reducing it's concentration to the ambient conditions.As stated earlier for the diffusive case (Section 3.2.2), the closed boundary is achieved byplacing an additional numerical cell beyond the model's end and equating the concentration in thiscell to that in the adjacent cell.Time Step and Grid Size Considerations As stated earlier, Ax and At will be given precedence over Cr, y, and Pe when deciding whichvariables to hold constant. C r , 7, and Pe are allowed to vary within the confines shown inFig. 3.6.Both Ax and At were chosen to be as small as possible to allow high accuracy and be capable ofmodelling the physical phenomena at hand. A second consideration is to chose Ax and At large38Chapter 3 - Description of Numerical Modelling39enough so that the total number of grid points and time steps do not exceed the computingcapabilities available.A time step of 1.25 s and a grid size of 5 cm (being the same order of magnitude as the channeldepth) were considered sufficiently small increments to model transport and mixing processesoccurring within the 1-D channel, while maximizing the computing power available. The 7.4 mchannel is therefore modelled using 148 grid points along its length and the typical 20 minute tidalcycle is modelled using 960 time steps.The final numerical model has been applied to the test cases studied within the physicalexperiments and presented in Table 4.1. The numerical results were acquired after the governingparameters were adjusted to provide the closest fit when compared to the physical case. Thesegoverning parameters are as follows:D = variable m2/s, Ax = 0.05 m, At = 1.25 s, Cr = variable, Pe = variable, andy = variableThe concentration results (both physical and numerical) are shown in Fig. 5.1 for the fiveconcentration plots ( xc = 1, 2, 3, 4, and 5) relating to the 'Base Case', over two tidal cycles.Truncation ErrorWhereas the truncation errors described for the previous cases contained only one or twovariables for the first, second, and third order-terms in E, the Taylor Series expansion for [3.33]becomes extremely complicated due to the presence of the area term, A r(x,t), and the varyingvelocity, U(x,t) in [2.9]. The associated truncation error contains terms of zeroth-, first- andsecond-orders implying that conservation of mass, advection, and diffusion will all be affected tosome degree within this FDM scheme. The important question is, how large are these errors andwill they significantly affect the results? The general form of E is:Chapter 3 - Description of Numerical Modelling40 t a2c Ax2 a3C E ate— n 1 + (UA - DAx) 31 ax3 + (higher order terms).4.: [3.37]Upon expansion of a 2ciat2 into a second-order spatial derivative, [3.37] is shown to contain 43variables in the zeroth order term, 29 variables in the first-order term, 11 variables in the second-order term plus other higher order terms.A sensitivity analysis was conducted for the grid size and time step used to determine the effectsof the zeroth, first, and second-order terms on the final solution. Although conservation of mass isinfluenced by zeroth order terms, the effects of the 43 variables in this situation are minimal. Thetotal mass of the pollutant within channel was altered by less than 10 -5% during the first two tidalcycles of testing.The first-order term was studied to determine its influence on the advective processes. Again,the 29 variables within this term only caused minor variations in the channel velocity (alteringU(x,t) by less than 10-4 % at any time). The second-order terms which influence the diffusiveprocesses were found to alter the diffusion coefficient by less than 0.01%.The reason that effects of such large groups of variables are so small is due to the fact that eachvariable is comprised of several high order derivatives of A r(x,t) and U(x,t) with respect to timeand space. These derivatives have orders of magnitude ranging from 10 -3 to 10-7 and whencombined, produce even smaller magnitude variables. The effects of E on C (or M), U, and D aretherefore not considered appreciable, but it would be unwise to believe that the solutions providedare without error.Numerical Tests The numerical tests run for this model used the same tidal conditions as used for the physicalmodels which are described in the following chapter. Input data for the hydrodynamic componentconsisted of tidal amplitude, A, channel length, X, average channel depth, d, and the tidalChapter 3 - Description of Numerical Modellingperiod, T. Figure 3.9 shows the velocity variations along the channel length during various stagesof the tidal cycle as depicted by the closed-form solutions for the Base Case. Input for theadvection/diffusion portion of the model consisted of the pollutant input position, xp, the initialconcentration C(x,0), and the initial discharge time, to/T.Figure 3.10 provides a 3-D perspective plot showing the path of the pollutant cloud in both timeand space. One can see how the pollutant is transported by both advection and diffusion. Thecross-sectional planes noted in this figure relate to positions within the physical models wherepollutant concentrations were measured. Figure 3.11 is an elevation view of the same datapresented in Fig. 3.10. The final data comparison is presented in the form of these 2-Dconcentration plots.41Chapter 3 - Description of Numerical ModellingChapter 4 - Description of Physical ExperimentsThe physical models for this research were used primarily to provide comparisons withnumerical results by measuring and correlating values of concentration, C(x,t) and velocity, U(x,t)for the desired channel configurations. These models are then used to refine specific diffusioncoefficients for each test scenario. All experiments were performed in the Hydraulics Laboratoryof the Civil Engineering Department at the University of British Columbia.Data acquisition and analysis was conducted using various computer platforms and videoequipment. Velocity and concentration measurements were made at locations corresponding topoints in the numerical grids of the computer programs. The number of simultaneous samplingpoints for a given test depended on the degree of accuracy needed and on time constraints. Sinceone of the limiting factors in the experimental set-up was the availability of only a singleconductivity probe, each test had to be repeated several times in order to ascertain what wasoccurring during a given test scenario.4.1 Testing ProcedureFigure 4.1 shows sketches of the channel used in the experiments. A pump connected to thetidal basin was used to control the channel's tidal rise and fall (flood and ebb). Tests wereconducted by fluctuating the water level at the basin end of the channel/basin structure. Once thetidal flow within the channel was established, effluent was injected at location x p and theconcentration and velocity over time were recorded for a specific point x c . The same test was thenrepeated measuring the concentration and velocity at a different position, x c . The average testscenario consisted of five individual tests referred to as test segments. A complete group of these42Chapter 4 - Description of Physical Experiments43tests segments is referred to as a test series. As shown in Table 4.1 the testing procedure for thephysical experiments consisted of varying the following parameters:1) Pollutant discharge position, xp2) Average water depth, d3) Tidal period, T4) Tidal amplitude, A5) Pollutant discharge position in tidal cycle, toaA Base Case was chosen around which these five parameters were varied. This base case is asfollows:xp = 3.5 m,d = 8.0 cm,T = 20 min,A = 30 mm,to/T = 0, 1, 2...(Peak Ebb) and 1/2, 3/2, 5/2,...(Peak Flood)During the modelling process all parameters were held constant except one. The significance ofthis lone parameter to the pollutant transport was then studied. The extent to which the parameterswere varied depended on limitations of the testing apparatus (as discussed herein).These parameters may be related directly to the dimensionless parameters developed in Section2.4. The observations from these tests will be discussed to a greater degree in Chapter 5 where theexperimental results are compared to those from the numerical model.4.2 Development of ApparatusWithin this section the various parts of the experimental set-up are described and relatedproblems discussed. The instrument set-up is shown in Figure 4.2.Test ChannelAs shown in Fig. 4.1, the one-dimensional test case was modelled in a narrow, constant width,flume which is 7.35 m long, 0.25 m wide, and 0.48 m deep, with a large rectangular basin whichChapter 4 - Description of Physical Experiments44is 1.07 m wide, 0.91 m long, and 0.49 m deep attached at one extreme to simulate the channelmouth (x = 0). The floor of the tidal basin was lowered 16 cm below the floor of the channel inorder to model the situation of deeper conditions offshore.Tidal Pump The pump configuration consisted of a variable speed peristaltic pump controlled by an inputvoltage from a signal generator which produced sinusoidal voltages of the necessary amplitude andfrequency. The maximum pumping capacity (13/ / min.) provided an average drawdown rate of 3mm /min. which determined the upper velocity limit for the test series. For example, the 20 minuteperiod chosen for the base case provides a maximum tidal amplitude of 30 mm and a subsequentmaximum based on A and T.During the test procedure pumping rates were found to fluctuate significantly due to thedeterioration of the flexible pump tubing. The abrasiveness of the pump head on the tube wouldpartially collapse the tube walls during testing, thereby reducing flow rates. Furthermore, smallpieces of abraded tube (made of a synthetic gum rubber) dropped into the pump head, slowingdown the rotating mechanism.As a result, regular calibration was necessary to determine the actual pumping rate for a givenoutput signal (delivered from the signal generator). These fluctuating flow rates complicated thetask of reproducing each test segment within a series five times. The tidal water flowing to andfrom the tidal pump passed through a series of diffuser tubes to ensure a uniform velocity profile inthe flume.The original pump diffuser apparatus - a single 12.5 mm dia. tube - caused jetting and undesiredcurrent patterns within the basin and channel. To alleviate this problem, flow was channeled - viaeight 5 mm tubes - to two horizontal rows of plexiglass tubing placed one atop the other andapproximately 10 cm apart. Each of the tubes was perforated with 3 mm dia. holes every 20 mmChapter 4 - Description of Physical Experiments45along its length. The two tubes were placed deep enough within the basin to avoid surfaceshearing but high enough to by-pass the deeper layer responsible for the upwelling at the channel'slip. To further ensure a nonstratified flow, the diffuser apparatus was placed behind a permeablemembrane (or buffer) made of perforated plywood covered with horsehair matting.Effluent Outflow DeviceThe effluent was diffused through a series of outflow tubes to minimize interference with thechannel's tidal flow. Momentum effects of the inflowing liquid from a single tube interrupted flowpatterns before the device was redesigned. The multiple tube array that was developed alsoallowed for faster effluent injection - in the order of 1.81 / min. of solution. The outflow tubeswere arranged in a grid fashion to simulate the inflow into a cell such as was modelled by thenumerical programs.This array consisted of four tubes extending into the water column forming two rows in line withthe channel and spaced equidistant from each other and the channel walls. Each tube had holesdrilled every 2 cm along its length, oriented at 90° to the flow direction. The device allowed thelongitudinal distance between tubes to be altered to allow for various input areas.Wave ProbeThe water surface elevation, 11 along the channel was measured using a capacitance-type waveprobe. This probe was designed by the National Research Council of Canada (N.R.C.) and builtat the University of British Columbia (UBC). The probe was positioned at the mouth of thechannel (x = 0) and became an essential monitor for adjusting flowrates from the somewhatunpredictable tidal pump. A calibration of the probe was conducted only once at the outset ofexperimentation, but visual checks were performed daily to ensure recorded results correspondedto channel depth markings.Chapter 4 - Description of Physical Experiments46The initial calibration showed that a linear relationship existed between the wave probe's outputvoltage and the water level. A calibration factor was therefore incorporated directly into theacquisition system and results were logged as water depth in mm from the channel's base. Resultswere accurate to the nearest 0.1 mm.Current Measurements Velocity measurements were obtained by visually noting the progression of dye tracers placed inthe channel at the same locations used by the conductivity probes. These traces consisted ofinjecting small amounts of rhodamine-wfrm dye into the channel with a syringe at mid depth.Measurements were taken every 60 seconds and the results were assumed to be representative forthat time period. Once more, measurements were taken over the central portions of the channel,outside the zone of influence of the boundary layers. This process was somewhat labour intensiveand carried out for certain test series (such as the Base Case) in order to allow a comparison withthe numerical model.Conductivity Probe During the testing period, three different probe layouts were used in conjunction with the threedifferent input positions, xp tested. These are shown in Figure 4.3. From the voltages measured,the corresponding salt concentrations were then determined. The probe works by measuring theconductivity of the solution across a 25 micron gap between two sets of platinum wires. Theconductivity probe corresponds to Type All manufactured by Precision Instruments (Head, 1983).The four platinum wires (dia. 5 ,um, each) at the tip of the probe are connected together in awheatstone bridge to provide the output voltage signal. Within these tests, the accuracy of theoutput data was dependent on the data acquisition system (DAS) and was set at 0.0025 volts. Theminimum discernible salt concentration associated with conductivity measurements within theseexperiments was 0.01 ppt.Chapter 4 - Description of Physical Experiments47The corresponding relationship between salinity and conductivity is indicated in Fig. 4.4 andmay be completely described by a polynomial (Lewis, 1980). However, the relationship isapproximately linear and may be described by a series of straight lines with slight changes in slope.The accuracy of this fit appears to increase as salinity approaches zero. This was beneficial sincemost salinity readings were less than 0.1 ppt.The output voltage was calibrated daily using a set of standard salinity solutions of the followingconcentrations: 0.0, 0.1, 0.2, 0.5, 1.0, 2.0, 4.0, and 10.0 ppt. It was noted that over a period of5 or 6 hours the probe's output voltage for a given salinity would increase. This "creep" isattributed to temperature variations caused by the instrumentation warming up. Whenevernoticeable changes occurred the probe was recalibrated.The conductivity probe was mounted at centre height within the channel to obtain a vertically -averaged concentration. The channel's central portion was least affected by boundary or shearlayers.To prevent damage to the extremely fragile probe tip, a tubular plastic housing was left in placeduring the experiments. This housing featured vertical slots extending the length of the contactzone of the probe. Nevertheless, the effect of this housing was to slow the flow field around theprobe's tip resulting in potentially inaccurate values.Data Acquisition SystemData acquisition was performed using a 12 bit Analog-to-Digital conversion board from KeithleyMetrabyte Corporation which was installed in a IBM compatible PC. The data aquisition programwas written in Quick Basic and made use of two of the eight available channels to record waterlevel data (calibrated within the program), and conductivity voltages (calibrated at a later date),along with the time of each observation. Recording intervals were set at one second. The programChapter 4 - Description of Physical Experiments48stored information in binary format which was subsequently converted to ASCII format foranalysis.Data sets from each test segment were later combined to produce data sets for each respective testseries using Lotus 123 speadsheet software. These experimental results were graphed andcompared visually with the numerical model. This is discussed in Chapter 6.4.2.1 EffluentProblems were encountered in producing a suitable neutrally-buoyant solution (specific gravityof 1.0) that could be traced using a conductivity probe. Following similar experiments describedby Fischer (1981), an effluent mixture of saltwater and methanol was considered. Common housesalt (NaC1, s.g. 2.17) breaks down readily in both water and methanol and provides an electricallyconductive solution. The methanol (s.g. 0.79) provides a density compensation for the heaviersalt. The addition of methanol poses few extra concerns since its vapour pressure is sufficientlylow to neglect short-term evaporation and quantities are small enough not to pose any health risk.Initial effluent mixtures used standard volumes and standard densities to produce a neutrally-buoyant solution. This mixture was not successful and quickly sank to the channel floor where itmoved as a density flow. After some time it was discovered that the solution should be mixedusing the molar volumes of the two liquids for the following reasons. When methanol is mixedwith water, a small percentage of the methanol dissolves in the water thereby leaving a smallervolume than the two original liquids (but with the original mass). By using molar volumes thisproblem was mitigated.Temperature Effects Although it was known that a neutrally-buoyant solution could be calculated, some densityproblems still persisted. Initially, temperature variations between the effluent (normally at roomtemperature) and the channel water (some 5 degrees cooler than room temperature) caused theChapter 4 - Description of Physical Experiments49solution to float. The solution was then kept in a cold-water bath (being the same temperature asthe channel water) to equilibrate the temperature of the two liquids. During the test program,temperature effects of this nature caused several tests to be aborted.Furthermore, the density effects of strong salt concentrations were more difficult to neutralizethan those of more dilute ones. It was, therefore, decided to use a relatively dilute solutioncontaining 10.0 parts per thousand (ppt), NaCl. Upon entering the channel, the effluents'concentrations dropped immediately to a third of the original value and rapidly decreased toapproximately 1.0 ppt, NaCl. Average long-term concentrations were in the order of 0.1 ppt,NaCl or one percent of the original concentration.Rhodamine-wt Dye To track the effluent visually during the testing program, rhodamine-wt dye was added to thesaltwater/methanol solution in small amounts. Because the amounts of dye added were extremelysmall, no density adjustments were made even though the specific gravity is equal to 1.14.Generally, solutions were mixed in the laboratory and their buoyancies were then tested in thechannel water by visual inspection. If a solution was not neutrally-buoyant (for whatever reasons)small amounts of methanol or salt solution were then added to compensate for inconsistencies.4.2.2 Aspects of FlowTest Initialization At the initiation of each test, the channel water was first cleared of residual currents by placingdividers within the channel at regular intervals. Dye droplets were added to the channel water andobserved to see if movement occurred. The test series then commenced with the tidal pumpproceeding through half a tidal cycle before the effluent was released. Both of these precautionswere taken to ensure that all currents were due solely to the tidal pumping and that steady stateconditions existed. A series of video images showing the movement of the pollutant cloud areChapter 4 - Description of Physical Experiments50presented in Fig. 4.5. Once the effluent had entered the channel, little else had to be done exceptfor measuring flow velocities and monitoring the water levels (controlled by the pump) to ensurethat they were indeed correct. If not, minor adjustments were then made to the pumping rates tocompensate.Boundary LayersA persistent difficulty with these experiments related to establishing a vertically uniform flow.Surface, wall, and bottom friction, in conjunction with the effect of the lip at the channel's mouth,and the pump's outflow configuration, all attributed to these vertical flow variations. A strongsurface shear layer was reduced using liquid dish soap to reduce surface tension effects. Wall andbottom friction were kept to a minimum by scrubbing the tank's sides daily to prevent a build-up ofsurface scum caused by the semi-stagnant channel water and the combination of substances withinthe water. Typical thicknesses of the boundary layers were approximately 10 mm at the surfaceand 15 - 20 mm along the bottom and sides and were fairly constant throughout the range of testsperformed.The 8 cm lip at the channel / basin interface caused upwelling and a subsequent surface shear onflood tides. This effect was considered an unavoidable feature of the experiments. Initially,horsehair matting was placed at the entrance in an effort to prevent the shear effects but this did notalleviate the problem. Wire mesh was placed across the channel at regular intervals which helpedreduce vertical variations in velocity and to induce some turbulence into the flow.A problem related to low flow velocities was that the corrresponding Reynolds number, Re (=Ud/v). This number, which is indicative of the extent of turbulent flow, is a function of depth,velocity and viscosity. A reynolds number of greater than Re = 10 5 is needed for a turbulentchannel flow (giving us a turbulent diffusion coefficient) (Henderson, 1966). During theexperiments however, the Reynolds number was generally less than 102 .Chapter 4 - Description of Physical ExperimentsChapter 5 - Results and Discussion5.1 General Observations from the Physical ExperimentsRepeatability of Results Since only one conductivity probe was available, each test was repeated five times, with theconductivity probe placed at different locations, x c , to ascertain the concentration levels occurringthroughout the channel. It was important to verify that conditions were similar during eachrepetition of the experiment and therefore several tests were repeated with the conductivity probe atthe same location. As shown in Fig. 5.1 (corresponding to the Base Case), measurements wereconducted twice with the conductivity probe at the same location, and the variations ofconcentration with time were noted. Although some variation does exist between the two curves,their trends and phases are reasonably similar. It is therefore assumed that a test series could berepeated with a satisfactory level of confidence and that conditions over the five repetitions did notchange greatly.5.1.1 Influence of Effluent Discharge Position, xpntxik = 0.20, 0.47 (Base Case) and, 0.73 During testing at Xpa = 0.20, significant portions of the effluent were quickly carried into thetidal basin during the following ebb tide where the effluent was subject to mixing and concentrationlevels on the following flood tides were too low to measure. This may be likened to the case oflongshore drift. This location also provided the fastest initial velocities to the newly-releasedeffluent and the shear layers in the effluent cloud for this case were less defined than in other tests.51Chapter 5 - Results and Discussion52For xp/X = 0.47 (i.e, the Base Case) the entire set of tests the transport of effluent was "wellbehaved" and no significant amounts of pollutant exited the channel within the two tidal cycles.After two tidal cycles the progression of the cloud in the seaward direction was noticeably furtherthan towards the landward end. Ebb and flood tides were transported equally well.Tests at XpA, = 0.73 (near the landward end of the channel) produced results that were much lessvaried than in the previous two cases since the effluent moved along the channel at a slower rate.At the slower velocities the wire mesh did not appear to aid in mixing and the dye cloud showeddistinct shear layers.5.1.2 Influence of Depth, d/kda = 0.01 and 0.02 Most of the tests were conducted at dA = 0.01 (Base Case) since it was considered that the depthto length ratio should be maximized in order to approximate actual estuarine proportions. This wasthe smallest depth attainable without having excessive bottom and wall friction. This ratiotranslates well into a typical 1-D channel and allows a greater A/d ratio to be tested. Havingd/A. = 0.02 modelled a situation with half of the effective length of the previously mentioned testThe associated reduction in tidal amplitude meant that the effluent transport was greatly reduced.Greater depths lead to less prominent shearing and a more consistent central layer since theboundary layers thickness stayed roughly constant.5.1.3 Influence of Tidal Amplitude, A/dA/d = 0.09, 0.19, and 0.38 (Base Case) Tests at A/d = 0.09 provided very little currents and the extent of transport after two cycles wasapproximately 2 m in either direction for the base conditions. The test series at A/d = 0.19produced stronger currents and a considerably more stable effluent case. The extent of effluenttransport after 2 tidal cycles was approximately 3 m. Due to the tidal pumping system available,Chapter 5 - Results and Discussion53Aid = 0.38 (equalling a 3.0 cm tidal amplitude) was the largest A/d ratio possible for a 20 minutetidal period. Transport after 2 tidal cycles reached the channel's mouth (x = 0 m). 5.1.4 Influence of Tidal Period, T2/g2,T = 20 min (Base Case) and 40 minAlong with the other factors involved, the 20 minute tidal period provided a situation where -after 2 tidal cycles - effluent was on the verge of exiting the channel. The 20 minute period wasalso the longest period appropriate with respect to the overall test schedule. At T = 40 min tests(for an amplitude of 3.0 cm) produced flows that contained large or predominant shear layers.Tests were too long to be considered practical. 5.1.5 Influence of Discharge Time, to/TEbb, to/T = 1.0 and Flood, taT = 0.5 In general, for to/T = 1.0 (discharging during the ebb), this situation allowed more of the effluentto enter the faster currents near the channel mouth. Although neither the flood or ebb tides exitedthe channel in general, mixing appeared to be more complete in the ebb case. For t o/T = 0.5,(discharging during the flood) the transport of effluent upchannel on the initial stage of inputcaused exit time within the channel to increase slightly.5.2 Combined Results of the Physical and Numerical ExperimentsAs described in Chapter 4 and shown in Table 4.1, the physical experiments consisted of eightseries of tests, each simulating various tidal conditions. These tidal conditions were subsequentlymodelled numerically as described in Chapter 3 and the two models were then compared on thebasis of concentration (and to a lesser extent velocity) information at the five conductivity stationsdescribed in Chapter 4. For the most part, comparisons were based on superimposingChapter 5 - Results and Discussion54concentration time histories of the measured results with several numerical results representingdifferent diffusion coefficients, D.Concentration ResultsGeneral observations of the two data sets suggests that the physical traces resemble numericaltraces with diffusion coefficients in the range of D = 5.0 x 10 -4 to 1.0 x 10-2 m2/s. Thesediffusion rates are approximately one to two orders of magnitude larger than that expected on thebasis of molecular diffusion suggesting that advective processes were involved. For the numericalmodel, several traces are shown in Fig. 5.2 to 5.7 relating to four different diffusion coefficients,namely:D = 0.5, 1.0, 2.0, and 5.0 x 10-3 m2/sEach physical test series was compared to this set of numerical traces and the appropriate numericaltrace was then chosen based on the criteria described below. These results are shown for five ofthe eight test series to determine the diffusion coefficient which produced the best fit for the varioustidal conditions. Comparison of the physical and numerical data sets was conducted visually onthe basis of the following criteria:1) similarity of initial concentration peaks at a given station (in amplitude, phase, andduration)2) coincidence (amplitude and duration) of flatter concentration trace sections whichnormally follow the initial peakOf the five physical traces (i.e. xc positions) for each test series, emphasis was placed on thecentral three traces as resemblances between the two data sets were strongest near the point ofrelease, Xp. Moreover, comparison was based on the concentration traces over the first of the twotidal cycles used since similarities at any station diminished with time. During later tidal cycles,inconsistencies in the physical data became more pronounced (see Chapter 4).Chapter 5 - Results and Discussion55The results of the comparison are shown in Table 5.1 which shows physical results matched to a`best fit' diffusion coefficient. The quality of the fit is also given (rated as Good, Fair, or Poor).As an example, the comparison for Test Series 2a is described below.Test Series 2aIn Figure 5.2, a diffusion coefficient, D = 2.0 x 10 -3 m2/s was chosen overall as providing thebest fit between the numerical and the physical traces. Starting with the uppermost graph,xc = 4.0 m, one can see that the initial peaks (t/T = 0.75) have very similar amplitudes but theduration of the physical peak extends approximately 0.4 tidal cycles, past that of the numericalpeak. This suggests that the pollutant cloud was not transported away from this position duringthe test but remained in the vicinity of the conductivity probe. A set of secondary peaks, located atapproximately UT = 1.75, shows a somewhat slower increase in the physical concentration butsimilar amplitudes and phases as the numerical data.At xc = 3.5 m, the initial peak (t/T = 0) in the physical data has a somewhat lower amplitude thanthe numerical value due to smoothing of the physical data. The phase of the physical peak is offsetdue to a slow pollutant release. The secondary peaks are of comparable amplitude but differ inphase by approximately 0.25 of a tidal cycle. The physical concentration trace tapers to zero afterone tidal cycle, whereas all numerical traces suggest concentration levels, C/Co to be in the rangeof 0.1 to 0.3. This suggests that the pollutant cloud had been transported away from the centralposition and did not return.At xc = 3.0 m, two physical traces are available showing the variability between physical testresults. Although curves #1 and #2 both display similar amplitudes (in both their primary andsecondary peaks), curve #1 appears slightly before #2 and declines at a slower rate suggesting alonger residence time at this position. In both cases, the numerical models underestimate theduration of the primary peaks. For the secondary peaks, both physical curves exhibit the sameChapter 5 - Results and Discussion56phase as the numerical curves but their amplitudes were less than those predicted by the numericalcurves.The quality of the fit was considered to be good despite the various shortcomings alreadymentioned. Variations in the physical results are considered to be due in part to the presence ofboundary layers in the physical experiments. Furthermore, the fact that each test series iscomprised of data from five separate test runs reduces the dependability of these results. Theinability of the numerical model to estimate the relative duration of the physical peaks might be dueto the plastic housing protecting the probe. This housing could potentially trap fluid during periodsof lower channel velocities, thereby overestimating the residence time at that location. From theresults, one cannot state with confidence that the pollutant cloud will be situated in a particularposition at a particular time. The combination of vertical variations in velocity and verticalvariations in dye concentration combine to transport the cloud either faster or slower thanpredicted.One can, however, follow the progression of the primary peak downchannel in Fig. 5.2 fromxc = 3.5 m to xc = 3.0 m. At this point the tides switch from an ebb to a flood, and the pollutantcloud moves upchannel towards xc = 3.5 m and xc = 4.0 m. For the concentration traces, onemay join the primary peaks of the graphs to provide a clearer picture of how the pollutant cloud ismoving along the channel with time.A similar analysis was performed for each of the test series and the results shown in Table 5.1.Generally, initial peak amplitudes for the physical experiments at all stations were within 20% ofthe input concentration suggested by the numerical results. In two of the six cases, secondaryphysical peaks at the discharge position, xp = 3.5 m, occurred well before that predicted by thenumerical models suggesting that the pollutant cloud is not transported as far as anticipated due to aslower advective process.Chapter 5 - Results and Discussion57In several cases concentration peaks arise in the physical results where none are predicted.Conversely, in some instances physical results taper off to zero concentration, where the numericalresults suggest that the pollutant cloud should be present.Velocity ResultsVelocity data consisting of results from measured dye traces and from the hydrodynamicequations were obtained for Base Case tidal conditions and are presented in a similar manner to theconcentration data. As stated previously, physical velocity data consist of noncontinuous pointmeasurements, averaged both in time and in the vertical profile of the test channel. Figure 5.8shows the corresponding velocity information at five locations (xc = 2.0, 3.0, 3.5, 4.0,and 5.0 m). Since the numerical velocity data are described by the hydrodynamic equation [2.37],using only the known parameters, no best fit analysis is required.As shown in Fig 5.8, there is little resemblance between the two sets of data at any givenposition over the two tidal cycles. Although the respective amplitudes are of similar orders ofmagnitude, there are large discrepancies in phase and duration. Moreover, measurements of thedye tracers tended to produce sporadic results with velocities oscillating between positive andnegative for each subsequent velocity reading (i.e. every 1/8 tidal cycle).These inconsistencies are associated with variations in the vertical velocity profile generallyattributed to boundary layers. The physical data obtained was therefore dependent on the depth atwhich the dye trace was observed. These results imply that the relative amount of advectionwithin the two models is different and therefore would lead to variations in the transport rates asshown in the concentration data. Overall, one can see that velocities decrease in magnitude as onemoves upchannel from the open channel boundary (i.e., from x c = 2.0 to 4.0 m).Chapter 5 - Results and Discussion58Overall ResultsFor the majority of the figures, a diffusion coefficient of 2.0 x 10 -3 m2/s was chosen. Thisrepresents a value greater than molecular diffusion (by approximately two orders of magnitude)suggesting that pollutant transport is therefore a function of both advective and diffusive processes.Within the physical model the shear layers present would create larger interfaces over whichmolecular diffusion occurs. Generally, one would expect that the diffusion coefficient should berelated to the velocities present within the channel and for tidal conditions which portray similarchannel velocities, the diffusion coefficient should be the same. Within the testing program, onlythree different velocity conditions were tested with five of those being similar to the velocityconditions for the base case. Following this line of logic, these test series should all match to onediffusion coefficient. As seen in Table 5.1, there is no perfect agreement with this and for testseries with similar velocity conditions as the Base Case (U max = 0.73 cm/s), D ranges from 0.5 to2.0 x 10-3 m2/s.As one can see, there is no set pattern relating maximum velocity to diffusion coefficient asanticipated. The variability in the physical results causes one to question the likelihood of makingan accurate comparison.5.3 DiscussionUsing a diffusion coefficient of 2 x 10 -3 m2/s concentration contours are presented in Figures5.9 to 5.12 showing the relative effect of the various tidal conditions to pollutant transport.Results are shown in terms of three of the dimensionless parameters developed in Section 2.4,namely:1) Tidal amplitude to depth ratio, A/d2) Initial pollutant discharge position, xp/A,3) Discharge time within the tidal cycle, tp/TChapter 5 - Results and Discussion59Tidal amplitude to depth ratio, A/dIn Fig 5.9, the tidal amplitude to depth ratio is varied between 0.05 and 0.19. This relates to25%, 50%, and 100% of the A/d ratio corresponding to the Base Case tidal conditions. One cansee that the results vary mostly in the area affected with larger A/d ratios and that after 2 tidal cyclesthe width of the pollutant cloud is the same for all three cases (approximately 0.3 of the channel'slength, X).The affected area, however, varies from 0.30 X for Fig. 5.9 (a) to approximately 0.47 X forFig. 5.9 (c). Also the peak concentration at any time is identical for all three cases, again owing toidentical diffusion coefficients.Initial pollutant discharge position, xLaIn Fig. 5.10, the pollutant discharge position, xp/X is varied between 0.20 and 0.74,corresponding to the discharge of the pollutant in regions of differing channel velocities. As seenin the figure, the range of the pollutant cloud decreases for the case of discharge positions furtherupstream, corresponding to lower velocities.In the case of Fig. 5.10 (a) the pollutant cloud reaches the seaward boundary where some of thepollutant is transported out of the channel in question and does not return. The total amount ofpollutant is therefore diminished and the subsequent concentration contours taper off at a faster ratethan in Figs. 5.10 (b) and (c).Discharge time within the tidal cycle, to[In Fig.5.11, the pollutant discharge time is varied for a central location within the channel.Fig. 5.11 (a) corresponds to a discharge at the peak velocity during an ebb tide (i.e, U = -Umax)while Fig. 5.11 (b) relates to a discharge on the flood tide (U = U max )• In these two case thespatial extent to which the pollutant is transported up and down-channel is approximately equal.Chapter 5 - Results and DiscussionChapter 6 - Conclusions and Recommendations6.1 ConclusionsBoth physical and numerical model studies have been conducted to lend insight to the subject oftransport and subsequent diffusion of a neutrally buoyant effluent instantaneously released intoone-dimensional tidal bodies.Numerical model developmentSeveral numerical models relating to various conditions of advection, diffusion, andhydrodynamic processes were developed. Within the hydrodynamic component closed-formsolutions for the velocity and surface elevation are developed for this one-dimensional, tidallyforced, constant depth case. These are applied in conjunction with the finite difference method(FDM) for determining pollutant concentration within the advection/diffusion components of themodels. These models were studied to determine the effects of varying the primary modellingparameters: Ax, At, U, and D. Secondary modelling parameters (namely, Cr, y, and Pe) are alsodiscussed. Results are compared to exact solutions wherever possible and discrepancies areexplained in terms of truncation errors associated with the numerical schemes employed.For the cases of pure advection or diffusion, it was determined that the primary modellingparameters could be arranged to obtain results approaching exact solutions. For theadvection/diffusion case, however, it has not been possible to state which mix of the primarymodelling parameters would produce the best results.60Chapter 6 - Conclusions and RecommendationsOne of the difficulties with the numerical model stems from At and Ax being assigned constantvalues throughout. To overcome this, the most reasonable approach appears to be to keep the Atand Ax values as small as possible to allow the increased accuracy (due to the FDM's convergentnature) to counter the model's varying accuracy.Physical model developmentPhysical models were developed and implemented on the basis of a suitable selection of thedimensionless parameters governing the situation. For a series of test conditions, pollutantconcentration time histories were obtained at five locations, along with velocity measurements forselected tests. Velocity information was considered too erratic to be used as a basis for comparisonto numerical results. Shear layers were present in many test cases which caused non-uniformitiesin the vertical direction. This implies that concentration measurements were, to some degree, depthdependent. If account is taken of vertical variation in concentration and velocity results areconsistent with the numerical results. The testing procedure was considered successful overall.Comparison of numerical and physical models Dimensional analyses were conducted and within the numerical and physical models tested, thefollowing factors were varied: effluent input position within the channel, xp/2t; release time withinthe tidal cycle, t/T; initial water depth within the tidal cycle, dA; tidal amplitude to depth ratio, A/d;and duration of the tidal period, T. Concentration and velocity data from the two models werecorrelated to derive a diffusion coefficient for the numerical model developed. A diffusioncoefficient of D = 2.0 x10 -3 m2/s (for Ax = 0.05 m , At = 1.25 s and y and Cr both varying)provides the best match between physical and numerical data.The final numerical model was then re-run for the selected diffusion coefficients and the effectsof the dimensionless parameters were discussed. Results indicate that while tidal amplitude willdetermine the extent of transport within a channel, the amount of diffusion after a full tidal cycle is61Chapter 6 - Conclusions and Recommendationssimilar for all cases. The transport is, however, skewed down-channel towards the seawardboundary and the skewness increases with larger tidal amplitudes. Similarly, for different inputpositions, the primary factor is the decrease in transport rates further up-channel.Pollutant released near the open boundary is transported by the greater channel velocities present,thereby having the opportunity to exit the channel. Pollutant released at the landward boundarytends to remain stationary, and is acted upon only by diffusive processes. The variation in theproportion of advective and diffusive processes along the channel is, perhaps, the largest singlefactor affecting the transport of a pollutant within this one-dimensional situation.6.2 Recommendations for Further StudyNumerical ModellingWithin the numerical modelling section, improvement could be made to the process by adaptingthe closed-form hydrodynamic solutions for the one-dimensional case to a width-varying case. Thefinite difference schemes used could be improved by implementing a more complicated FDMscheme such as the GREATEST scheme developed by Leonard (1979). Further study of the roleof the primary modelling parameters within the advection/diffusion model would also be useful.Physical ExperimentsIn the case of physical model tests, the apparatus used was well designed although theperformance limitations were not anticipated. These studies could easily be continued with theaddition of the necessary equipment. The channel used was a manageable size and has thepossibility to provide further valid test results. A series of conductivity probes will allow theresearcher to carry out a complete test without needing to repeat an experiment several times.Conductivity probe measurements would be taken at various depths within the channel to62Chapter 6 - Conclusions and Recommendationsdetermine the effects of stratification. A more accurate method of measuring channel velocities isalso needed.The most crucial point to further physical studies would be the acquisition of a larger capacitypump. The scope of the study was severely impeded by lesser ideal flow rates and tidal amplitudesresulting in low Re values. Stratified flows may be minimized by setting the bottom of the tidalbasin flush with the channel bottom. Shear and boundary layers due to wall and bottom frictionwill continue to be a problem that will most easily be solved by using a larger testing tank.63Chapter 6 - Conclusions and RecommendationsReferencesAbbott, M. (1979) " Computational Hydraulics - Elements of the Theory of Free SurfaceFlows ". Pitman Press, Boston.Abbott, M. and Basco, T. (1989) " Computational Fluid Dynamics ". Longman Scientific &Technical, United Kingdom. Ch. 3 - Ch. 5.Bienfang, P. (1980) " Water Quality Characteristics of Honokohau Harbour: A SubtropicalEmbayment Affected by Groundwater Intrusion". Pacific Science, Vol. 34, pp.279-291.Chu, W., Barker, B., Akbar, A. (1988) " Modeling Tidal Transport in the Arabian Gulf ".J. Waterway, Port, Coastal and Ocean Eng. Div., A.S.C.E., Vol.114, pp. 455-471.Chu, W., and Yeh, W. (1985) " Calibration of a Two-Dimensional Hydrodynamic Model ".Coastal Engineering Vol. 9, The Nederlands, pp. 293-307.Fischer, H., List, E., Koh, R., Imberger, J., Brooks, N. (1979) " Mixing in Inland and CoastalWaters". Academic Press, San Diego.Fischer, H. (1981) " Transport Models for Inland and Coastal Waters ". Academic Press, SanDiego.Head, M.J. (1983) " The Use of Miniature Four-Electrode Conductivity Probes for HighResolution Measurement of Turbulent Density or Temperature Variations in Salt StratifiedFlows ". Ph.D. Thesis, University of California, San Diego.Henderson, F.M. (1966) " Open Channel Flow " . Macmillan Publishing, New York.Koutitas, C. (1970) " Mathematical Models In Coastal Engineering ". Pentech Press, London.Leonard, B.P. (1979) " A Stable and Accurate Convective Modelling Procedure Based onQuadratic Upstream Interpolation". Computer Methods in Applied Science and Engineering.Vol. 19 pp.59-98.64Lewis, E.L. (1980) " The Practical Salinity Scale and Its Antecedents" IEEE Journal of OceanEngineering. Vol. 5, pp.3-8.Nece, R. (1984) " Planform Effects on Tidal Flushing of Marinas". J. Waterway, Port, Coastaland Ocean Eng. Div., A.S.C.E., Vol. 104, pp. 251.Nece, R., Richey, E. P. (1972) " Flushing Characteristics of Small-Boat Marinas". Proceedings13th Coastal Engineering Conference. Vol. 3, pp. 2499-2513.Sarpkaya, T. and Isaacson, M. (1981). "Mechanics of Wave Forces on Off-shore Structures" .Von Nostrand Reinhold, N. Y., pp. 198-202.Schwartz, R.A., Imberger, J. (1988) " Flushing Behaviour of a Coastal Marina". Proceedings21st Coastal Engineering Conference., pp. 2626-2640.van de Kreeke, J. (1983) " Residence Time: Applications to Small Boat Basins ". J. Waterway,Port, Coastal and Ocean Eng. Div., A.S.C.E., Vol. 103, pp. 416-430.65Tables67Cr U (m/s) Ax (m) At (s) Profile #0.5 0.5 5 5 1, 2, 3, and 410 10 1, 2, 3, and 420 . 20 'v 1 2 3 and 4 0.75 .. 0.75 , 5 5 1, 2, 3, and 410 10 1, 2, 3, and 420 20 1, 2, 3, and 41.0 1.0 ' 5 5 1, 2, 3, and 410 10 1, 2, 3, and 420 20 1, 2,3,and4 1.5 1.5 5 5 , 1, 2, 3, and 4 10 10 1, 2, 3, and 420 20 1, 2, 3, and 4Table 3 .1 Numerical modelling testing program for a pure advection caseD (m2/0 Ax (m) At (s) Profile #0.10 0.10 5 1.25 1, 2, 3, and 410 5 1, 2, 3, and 420 20 1, 2, 3, and 40.17 0.17 5 1.25 1, 2, 3, and 410 5 1, 2, 3, and 420 20 1, 2, 3, and 40.25 0.25 5 1.25 1, 2, 3, and 410 5 1, 2, 3, and 420 20 1, 2, 3, and 40.45 0.45 5 1.25 1, 2, 3, and 410 5 1, 2, 3, and 420 20 1, 2, 3, and 40.50 0.50 5 1.25 1, 2, 3, and 410 5 1, 2, 3, and 420 20 1, 2, 3, and 40.55 0.55 5 1.25 1 2 3 and 410 5 1, 2, 3, and 4-rrwr.-n•-•-•,■-■■■■■-■-■-■-■-■,■+■-■■■■-r.l 20 1, 2, 3, and 420Boundariesopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closed open closedclosed ' closedopen closedclosed closedopen closedclosed ' closedopen closedclosed ' closedopen closedclosed ' closed68Table 3.2 Numerical modelling testing program for a pure diffusion case0.0840.167 1111111111 133.6 400 133.6 400 133.6 400 200 133.6 400 2000.4900.666 2.000.7450.8330.490SIM0.735 1.500.6130.8582.501.0069Cr Pe0.084 1.00.167 2.003.000.250 0.333 4.000.4180.167 1.000.334 2.000.500 3.000.584 3.500.668 4.00D (m2/s) Ax (m) At (s)67.2 400 20067.2 400 20067.267.2 4004001111111111M1111266.0 400 200266.0 400 200266.0 400 200392.0 400 200392.0 400 200392.0 400 200392.0 400 200392.0 400 2000.980 2.00Numerical modelling testing program for Advection/Diffusioncomponent with constant velocityTable 3.3TestSeriesXp,DischargePosition(m)d,AverageWaterDepth(cm)T,TidalPeriod(min)A,TidalAmplitude(cm)Umax,MaximumChannelVelocity(cm/s)DischargeTimeFloodMEMIN 20 .01 10.0OM MI0.38 0.18 Ebbb =III 111111111111111111111111111111 Flood4 a 1111111111111111 20.0 0.75 0.36 Ebbb 1111111111111111111111111.11111111111111 FloodIMIIIIIIMM 8.0 20.0 11 Ebbb 1111111111=1111111111111111111111111 Flood6 a 11111.1 40.0 3.0 0.73 Ebbam ffilliall111111111111b IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIM Flood8 ab)MEM 8.0 20.0 NMI 0.73 EbbFloodNote: The highlighted lines refer to the Base CaseTable 4.1 Physical testing programme showing variation in tidal conditionsbased upon xp, d, T, A, and the discharge time.7071Test Seriesxnr 'DischargePosition(m)Av gageWaterDepth(cm)'I', TidalPeriod(min)A ,TidalAmp.(cm)Umax,Max.Velocity(cm/s)DiffusionCoefficient,D(10 -3 m2 /s)Quality ofFit.2 a (Ebb)- --3.5 8.0 10.0 0.75 0.73 1.0 Goodb (Flood) 1.0 Good3 a KEbb) 3.5 8.0 20.0 0.38 0.18 2.0 Fairb (Flood)_ , 1.0 - 2.0 Fair4 a (Ebb) 3.5 8.0 20.0 0.75 0.36 0.5 - 1.0 • Poorb ood5 a (Ebb) 3.5 .... 8.0 20.0 1.5 0.731.00.5 - 1.0 ..,FairPoor b (Flood) 1.0 Fair .7 a bbb (Flood)3.5 16.0 20.0 : 1.5 0.36 2.0 - 5.0 2.0 - 5.0 .FairFair8 b (Flood) 5.5 8.0 20.0 1.5 0.73 2.0 FairTable 5.1 Best Fit Analysis comparing physical concentration traces with numericallyproduced traces for various diffusion coefficients.Figures72x(b) One-dimensional, variable width case-./ / / / / / / / / / / / / / /._,„„... x< ). < ). CLW) < y. j ..--(a) One-dimensional, constant width case(c) Two-dimensional case73Figure 1.1 Sketches of (a) 1-D, (b) 1.5-D, and (c) 2-D channelsx = 0 )1Effluent Inflowx = X74L(a) tidal notationAXv.".-SNN.N-NNNN.NvNvNNN'sN.V.NNNW.:.s.NNNNNNNN.NWNNNN\Ns.:s.NWsN'NNNNN'N*NNs.,b) elevation(c) planx = 0 Figure 2.1 Definition sketches of the one-dimensional situation(a) tidal notation (b) elevation (c) plan75Additional cell Additional cell I 1 \\\x = - Ax :::\r r\-\x\t\\\\\t\\\2\)N-\\\\\\\t\N\N\\ \ 1Numerical cell x = X-4-Ax I IOpen (seaward) boundary Closed (landward) boundary Figure 3.1 Layout of the finite difference scheme showing boundary conditionsand numerical cellsC x = 0—)-1 Ax k--- x = X1 point trianglet= 0Profile #11.Q76C/C,0.01 .0C/CoProfile # 23 point trianglet = 00.0 x )1, 1.00.0Profile # 3C/C,Profile # 43 point rectanglet = 0x1 .0C/Co0 . 0 .Note: the effluent for each model was equilibrated to contain a unit volume.Figure 3.2 Initial shapes of pollutant concentrations testedin numerical models.Exact Solution t= 0 t= X/UCr = 1.001.0 -0.5 -Cr = 0.50 - t = 0-- t= X/U770.00.00 0.25 0.50 0.75 1.00Cr = 0.75 - t = 0- - t= X/U1 .0]0.50.00.00 0.25 0.50 0.75 1.000.00 0.25 0.50 0.75 1.001.0 -0.5 - Cr = 1.500-t=t=0.05 ?/U t=0.10 X/U0.0-0.5 --1.0I I I0.00 0.25 0.50 0.75 1.00Length along channel, x/ANote: X/ U represents one tidal cycleFigure 3.3 Effect of varying the Courant number (= Udt1 Ax) within theadvection model.1.0 —0.5 —0.0Ax = 5 mAt =5 s t = 0--0.00 0.25 0.50 0.75 1.001.0 —Ax = 10 mAt = 10 s Exact Solution- t = 0- - t= X/U0.5 —C.)0.00.00 0.25 0.50 0.75 1.001.0 —Ax = 20 m -t= 0--At = 20 s0.5 —0.0 I I I0.00 0.25 0.50 0.75 1.00Length along channel, x/X,Note: X/U represents one tidal cycleFigure 3.4 Effect of varying At and Ax while holding the Courant number (= UAt/Dx)constant within the advection model.7879— t=Os t=200s----t= 400 st= 600 sAx = 20 mAt = 20 s50-5-101.0 —= 0.10 InitialPollutantShapeExact Solution,t = 600 s0.5 —• - -------------0.050-5 10-1050 10-10 -5-5 10-10 0 5y=0.17Exact Solution,t = 600 s0.0----------------0.0= 0.55Exact Solution,t = 600 s11-0.5 —1.0 —(-) 0.5 —U1.0 --0.5 —= 0.45Exact Solution,t = 600 s1.0 —0.5 —0.0Gridspaces from initial position, xpFigure 3.5 Effect of varying y(= DAtlAx2 ) within the diffusion model byvarying D while holding At and Ax constant1.00.9 -0.80.70.6Cr0.50.40.3 -100.20.11.00.50.2-J-0.1 0.2 0.3 0.4 0.5Figure 3.6 Stability range for FTCS finite difference operator with respect toCr and y(after Leonard, 1979)800.00 0.25 0.50 0.75 1.000.00 0.25 0.50 0.75 1.000.5 —0.0-0.5 —t = 0Pe = 1.0Exact Solution0.5 —0.01.0 —t = 0Exact SolutionPe = 2.01.0 -,(2 0.5-U0.01.0 —t = 0Exact SolutionPe = 4.00.00 0.25 0.50 0.75 1.00Length along channel, x/XFigure 3.7 Effect of varying Pe (= C In within the advectionldiffusion modeland varying C. Ax, and At while holding 7, U and D constant81S.Pe = 5.0/..Pe = 4.0.Pe = 2.0-= 1.0 ---0.75 1.000.50\ ..Pe = 2.0= 1.5--..----Pe =0.25I I . I0.50 0.75 1.00Distance along channel, )(A.= 0Pe = 3.5,-Pe =....Pe = 1.0....... ......0.50 0.75 1.00t = 0Figure 3.8 Effect of varying y(= DAtlAx 2) within the advection/diffusionmodel while holding dr and At constant1 .0 -0.5 -0.0 = 0.084S.S.0.00 0.251.0 -y = 0.167u° 0.5 -*C.D.... /0.0 -0.00 0.251.0 -S.'Y= 0.4900.5 -0.00.0082t.0g 1.00 -E0.67t =0, 1/2 T, T3/4 Tt = 5/8 T,7/8 Tt = 1/8 T, 3/8 Tt = 1/4 T79-0.00 0.25 0.50 0.75 1.00Distance along channel, x/XSeaward End Landward EndFigure 3.9 Velocity variations along channel length during various stages of the tidal cycleas predicted by closed-form solutions for the Base CaseNI NIN I NNNN1 .0NN NC/CoN0.0 N-24 gitiott.-4" ItziK44,NO, -0004,,etottea.`NN1.0`GPIesx =Pxc = 0.25(1,,, /19'2.0 e0x = zysa,xc = 0.74 A,x = 0.47 A.Cross-sectional planes through whichconcentration traces are takenNNNgyFigure 3.10 3-D perspective plot showing path of pollutant cloudin both time and spaceDischarge Position, xp0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 11.000.900.800.700.600.500.400.300.200.100.0085Distance along channel, xASeaward End Landward EndFigure 3.11 Concentration contours in time and space for the Base Case.(Discharge on the ebb tide)Concentration MonitorsServo-Controlled PumpInflow / Outflow Effluent InflowDiffuser Wall1-D Channel7.35 m0.48 mitxp TidalI — BasinHo- 0.91 mVolume Channel 0.90 m 3Volume Basin 0.48 m 386Figure 4.1 Plan and elevation of channel used in the physical model studies.EffluentReservoirConductivity Probes•TidalPumpTidal ReservoirPowerSupplyVoltage ConverterSignalGeneratorOscilloscopeProbePreAmpMicrocomputerDataAcquisitionSystem87r -Figure 4.2 Set-up of physical model instrumentationWave probe Conductivity ProbesXp-►dx = 1.5 mate =-0.0, 1.0,1.5, 2.0, and 2.5 m1-01- XP x = 3.5 mXc = 2.5, 3.0, 3.5, 4.0 and 4.5 mx_P = 5 5 mxc = 4.5, 5.0, 5.5, 6.0, and 6.5 mFigure 4.3 Position of conductivity probes, xc, in relation to variouseffluent input positions, xp8889-3-450I I I I I I1 2 3 4 5 6Salinity, ppt1 7 8 9 10I IA — Measured. ConductivityFigure 4.4 Correlation for the conductivity probe between measuredconductivity and the predetermined salinity values• 0.50 0.55............... - . : . • • :maw 'q o0.450.35 0.40 0.650.600.000.08.E 0.15g0.23OE 0.30z0.38Open (seaward)boundaryDistance along channel, x/X.Closed (landward)boundaryFigure 4.5 Video images showing the movement of the pollutant cloud forthe Base Case. (Input on the flood tide)1.00.8Physical Results0.6 Test Series 5 i(b)a) Test Series 5 (b) iiOU.g 0.40.20.000ISlack TideU = 0 0.5' Slack TideU = 0 1.0Slack TideU= 0 1.5Peak Ebb Peak Flood Peak Ebb Peak FloodU = U = Umax U = -Umax U = UmaxTidal cycles, t/TFigure 5.1 Repeatability of effluent concentrations measurementsfrom physical tests for the Base Case.910.0 0.5 1.0 1.5 2.0xc = 3.5 mz 1 .0 1.50.50.0 2.01.5 2.00 0 0.5 1.0Tidal cycles, t/T• 0.01.0= 3.0 m0.80.6Curve #1Curve #20.20.00.41.0 -0.8 -0.6 -0.4 -0.2 -xc = 4.0 mInitial peakSecondary peak0.0.....; - • _..... . - - "0.2Diffusion Coefficient, D (10 -3 m2/s) 0.5- - 2.0- - - - 5.0- Physical Results1 .0O 0.8-4• 0.68o 0.492Figure 5.2 Concentration traces at locations x, = 3.0, 3.5, and 4.0 m comparingphysical and numerical results corresponding to Test Series 2a.(Discharge at xp = 3.5 m on the ebb tide, to IT = 0.0)x, = 3.5 m12.00.5 1.0 1.51.0 -C.)o 0.8te! R0.6 -C)0 0.40.2-Ri •0.00.0Diffusion Coefficient, D (10 -3 1112M 0.51.0- - 2.05.0- Physical Results1.00.80.60.40.20.00.0 0.5 1.0 1.5 2.01.00.80.60.40.20.00 0 0.5 1.0 1.5 2.0Tidal cycles, t/T93Figure 5.3 Concentration traces at locations x e = 3.0, 3.5, and 4.0 m comparingphysical and numerical results corresponding to Test Series 2b.(Discharge at xp = 3.5 m on the flood tide, toff = 0.5)0.8 -0.6 -0.4 -0.2 -0.0 •z1.00.8 0.5.630.6a)00.4c.)0.20.0= 3.5 mDiffusion Coefficient, D (10 -3 rn2A)0.0 0.5 1.0 1.5 2.0- - 2.0- - - - 5.0- Physical Results1 .0 -94x = 4 0 mT. • • •... ........ • ...............0.0 0.5 1.0 1.5 2.01.00.80.60.40.20.00 0 0.5 1.0 1.5 2.0Tidal cycles, tITFigure 5.4 Concentration traces at locations x, = 3.0, 3.5, and 4.0 m comparingphysical and numerical results corresponding to Test Series 3a.(Discharge at xp = 3.5 m on the ebb tide, tolT = 0.0)1.00.80.60.40.20 .0950.0 0.5 1.0 1.5 2.0Diffusion Coefficient, D (10 -3m2/s) 0.5 - - 2.0--- 5.0- Physical Results ... I 1.5 2.01.00.80.60.40.20.00.0 0.5 1.0 1.5Tidal cycles, t/TFigure 5.5 Concentration traces at locations xc = 3.0, 3.5, 4.0 m comparingphysical and numerical results corresponding to Test Series 4a.(Discharge at xp = 3.5 m on the ebb tide, tolT = 0.0)2.0x = 3.5 m0.5 1.0 1.5 2.0C"r--5'' 1.0E 0.8 -C 0.680 0.4 --c10.2 -cd _0.0o.oDiffusion Coefficient, D (10 -3 rn2 Is) 0.5- - - - 1.0- - 2.0- • - 5.0- Physical Results1.00.80.60.40.20 .000 0.5 1.0 1.5 2.01.00.80.60.40.20.00 0 0.5 1.0 1.5 2.0Tidal cycles, tfT96Figure 5.6 Concentration traces at locations x, = 3 .0, 3.5 , and 4.0 m comparingphysical and numerical results corresponding to Test Series 5a.(Discharge at xp = 3.5 m on the ebb tide, tolT = 0.0)xc = 3.5 m1.5 2.00.5 1.0C.)) 1.00 0.8cc.t0.61.)800.2o 0.00.0Diffusion Coefficient, D (10 -3 1112/s) 0.5- - 2.0- Physical ResultsCurve #1Cue #240#1.4 1 11 11111111&11rv.-1.00.80.60.40.20.09700 0.5 1.0 1.5 2.01.00.80.60.40.20.00 0 0.5 1.0 1.5 2.0Tidal cycles, VI'Figure 5.7 Concentration traces at locations xc = 3.0, 3.5, 4.0 m comparingphysical and numerical results corresponding at Test Series 5b.(Discharge at x, = 3.5 m on the flood tide, t 0IT = 0.5)10.5A11.51.0 2.02.0. I . I I0 0 0.5 1.0 1.512.0AI ' 1 I0 0 0.5 1.0 1.500AA I I 0.5 1.0 1.5AA A A AAI 1 1 10.5 1.0 1.5 2.0A0.0A0.012.0-- A A---A, A- Ax = 5.0 m0.750.00-0.75X = 4.0 m0.750.00-0.75x c =3.5m ,,,--_-; 0.750.00. -8 -0.7571 )')x = 3.0 m 0.750.00-0.75xc = 2.0 m 0.750.00-0.75Tidal cycles, UTA - Measured VelocityFigure 5.8 Velocity profiles comparing physical and numerical resultsat various locations for the Base Case.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1020.511.52 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.51Distance along channel, x/A.Figure 5.9 Numerical concentration contours comparing various tidal amplitudeto depth ratios, Ald = (a) 0.05 , (b) 0.09, and (c) 0.199100.510 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 122 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Distance along channel, x/X,Figure 5.10 Numerical concentration contours comparing various input positions alongchannel length, VA, = (a) 0.20, (b) 0.47, and (c) 0.74.I co00.511.52 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.511.52Distance along channel xAFigure 5.11 Numerical concentration contours comparing various input times withinthe tidal cycle, tpIT = input on the (a) ebb tide and (b) flood tide.to I