MODELLING OF T H E CIRCULATION OF NORTHERN BRITISH COLUMBIA WATERS By Renee Jacques B. Sc. (Physics) Universite de Sherbrooke, 1991 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES E A R T H A N D O C E A N SCIENCES (OCEANOGRAPHY) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A May 1997 © Renee Jacques, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Earth and Ocean Sciences (Oceanography) The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: Abstract A numerical model is used to improve our understanding of the mean circulation of northern British Columbia waters. The velocity field arising from tidal, wind, river runoff, and baroclinic effects is studied and the influence of the earth's rotation investigated. The model is three-dimensional, diagnostic, nonlinear through tidal forcing, and uses finite element techniques. The grid has linear horizontal triangular elements with vertical sigma levels. The model output is a linear combination of the currents obtained under the individual forcings. In this thesis, the model has been modified in order to improve the parameterization of friction processes, to add the vertical advection of tidal momentum by the tides, and to include the stratification effect on eddy viscosity. A major circulation feature is a large cyclonic eddy covering the central-eastern Dixon Entrance. This eddy, referred as the Rose Spit Eddy, is actually simulated by a combination of two or three smaller gyres, which seasonally vary in size and intensity. The modeled eddy is, in good agreement with observations, characterized by a rotation period of 3 to 40 days, depending on the season and the initial distance of the numerical particle from the center of the gyres. This eddy persists in all seasons but occasionally disappears. A disappearance occurred, consistent with current observations, in October 1984. The Learmonth Bank Gyre, an anticyclonic eddy located in western Dixon Entrance, is also simulated by the model. Baroclinicity and tidal rectification, respectively responsible for the generation of the Rose Spit Eddy and Learmonth Bank Gyre, are the dominant effects in the study region. The earth's rotation, necessary for the formation of those two eddying motions, as well as the vertical advection of tidal momentum were found to be essential to bring the current ii speeds closer to the observed values. The effect of stratification on the residual currents, through turbulence reduction, was found to be negligible. 111 Table of Contents Abstract ii List of Tables x List of Figures xii Acknowledgements 1 2 xxiii Introduction 1 1.1 Motivation 3 1.1.1 Fisheries 3 1.1.2 Oil exploration 4 1.1.3 Pollution discharges 5 1.2 Objectives 5 1.3 Overview of Thesis 6 Review of the Physical Environment and Earlier Models 8 2.1 Topography and Bathymetry 8 2.1.1 Dixon Entrance and adjacent seaways 9 2.1.2 Chatham Sound 11 2.2 Sea Floor Sediments 12 2.3 Meteorology and Climatology 13 2.3.1 Winds 14 2.3.2 Precipitation 15 iv 2.4 2.5 3 Hydrology 15 2.4.1 16 Skeena and Nass rivers Physical Oceanography 17 2.5.1 Adjacent oceanic currents 18 2.5.2 Tides and tidal currents 19 2.5.3 Non-tidal currents 22 2.5.4 Gyral circulation 29 2.5.5 Water properties 30 2.6 Climatic Variability 32 2.7 Physical Oceanographic Models 33 2.7.1 Laboratory model 33 2.7.2 Barotropic numerical models 34 2.7.3 Baroclinic numerical models 36 Description of the Numerical Model and its Physics 38 3.1 38 3.2 Boundaries of the Model 3.1.1 Ocean floor 39 3.1.2 Free ocean surface 42 3.1.3 Land and open boundaries 42 Model Friction 43 3.2.1 Friction at the bottom 44 3.2.2 Friction at the surface 46 3.2.3 Friction in the water column 50 3.3 Effect of Stratification 54 3.4 Governing Equations of the Model 56 3.5 Numerical Solution Procedure 58 v 4 Solving the governing equations 58 3.5.2 Relaxation and convergence test 60 Modelling the Dixon Entrance System 62 4.1 Spatial Discretization 62 4.1.1 Horizontal discretization 63 4.1.2 Vertical discretization 66 4.2 Open Boundaries 71 4.3 Model Parameters 72 4.3.1 Friction parameters 72 4.3.2 Other parameters 74 4.4 4.5 4.6 5 3.5.1 Tides 75 4.4.1 Tidal constituents 76 4.4.2 Calculation of tidal constituents 77 4.4.3 Tidal forcing 78 Baroclinic Forcing 80 4.5.1 Density observations 80 4.5.2 Interpolation of density data 81 4.5.3 Calculation of baroclinic pressure gradients 87 Wind and River Discharge Forcings 93 4.6.1 Wind stress 93 4.6.2 River runoff 94 Tidal Rectification and Wind-Driven Circulation 100 5.1 Tidal Residual Currents 100 5.1.1 Generation mechanisms 101 5.1.2 The vorticity equation 102 vi 5.1.3 5.2 5.3 104 Role of the Vertical Advection of Tidal Momentum by the Tides Ill 5.2.1 The momentum equation Ill 5.2.2 The vorticity equation 112 5.2.3 Numerical results . 115 Wind-Driven Currents 117 5.3.1 Idealized wind 117 5.3.2 Observed wind 120 5.4 Influence of the Earth's Rotation 124 5.5 Sensitivity of the Model Results 127 5.5.1 Open boundary conditions 127 5.5.2 Grid resolution 129 5.5.3 Friction parameters 134 5.6 5.7 6 Numerical results and analysis of vorticity terms Stratification Effect 137 5.6.1 Modification to the governing equations 137 5.6.2 Numerical results 139 Summary 140 Currents Arising from River Runoff and Baroclinic Effects 143 6.1 Currents Driven by River Runoff 143 6.2 Currents Driven by Baroclinic Pressure Gradients 144 6.2.1 Summer currents 145 6.2.2 Fall currents 153 6.2.3 Spring currents 153 6.2.4 Winter currents 159 6.3 Model Sensitivity 159 vii 6.4 Influence of Earth's Rotation and of Stratification on Eddy Viscosity 7 6.5 Summary 161 Evaluation of Model Accuracy 164 7.1 Mean Water Circulation Driven by all Forcings 164 7.1.1 Numerical results 165 7.2 Comparison of Model Results with Current Measurements 8 . . 161 168 7.2.1 The 1962 current observations 168 7.2.2 The 1984-85 current observations 171 7.2.3 The 1991-92 current observations 176 7.3 Summary 181 Conclusion 182 8.1 Future Work 185 Bibliography 186 Appendices 197 A Glossary of Commonly Used Symbols 198 B Historical Review of Observations 202 B.l Sea Floor Sediments 202 B.2 Winds 203 B.3 River Discharge 206 B.4 Oceanographic Cruises 207 B.5 Current Measurements 209 B.6 Lighthouse Observations 211 Vlll C Governing Equations 212 C.l Continuity Equation for an Incompressible Fluid 212 C.2 Navier-Stokes Equations 213 C.3 Reynolds Equations 213 C.4 Reynolds Equations in the Earth's Rotating Frame 215 C.5 Boussinesq Fluid 216 C. 6 Shallow-Water Equations 216 D Derivation of the Model Equations 218 D. l Representation of Mixing Processes 218 D.2 Bottom and Surface Boundary Conditions 219 D.3 Vertically Integrated Equations 220 D.4 Frequency Domain 221 D.4.1 Absolute value term 222 D.4.2 Governing equations 223 D.5 Steady Component 224 D.6 Rotary Components 226 D.7 Partial Solutions to Momentum Equation 227 D.8 Momentum Equation Reformulation 228 D.9 Generalized Wave Continuity Equation 229 E Finite Element Method 231 F Statistical Interpolation 234 F.l 234 Basic Theory F.2 How to Apply the Theory to a Physical Problem ix 238 List of Tables 3.1 Neutral drag coefficients according to various authors (from Blake, 1991). 49 4.1 Model parameters 75 4.2 Tidal constituents, in order of importance, included in the model 76 4.3 Interval between two consecutive z-levels. The interval depends on position in water column, the water surface being at 0 m 4.4 81 Wind velocity and stress averaged over two cruise periods, June 24-July 7 and December 2-18 1991. The time average is given for the three ocean buoys shown on Fig. B.l as well as their spatial average. Wind direction is defined as the direction towards which the wind blows with angles measured clockwise from true north 4.5 98 Time-averaged discharge rate of the Skeena and Nass rivers computed over each of the cruise period 99 5.1 Description of the terms of the vorticity equation (5.5) 104 5.2 Description of the terms of the steady-state vorticity equation (5.6). . . . 110 5.3 Absolute value of each term of the steady-state vorticity equation (5.6) averaged over the five regions shown on Fig. 5.4. Units are 10 s . . . . 110 -lo 5.4 -2 Absolute value of term HH of the steady-state vorticity equation (5.15) averaged over the five regions shown on Fig. 5.4. Units are 10 s . The -10 -2 % represents the relative importance of term HH with respect to the most important vorticity term (bold value in Table 5.3) x 114 B.l Description of land wind station sites shown on Fig. B.l. At Rose Spit, an automatic wind sensing station is operated where wind direction and speed are radioed several times per hour to Prince Rupert. (From Atmospheric Environment service (1981; 1982) and Brower et al. (1977).) 205 B.2 Major oceanographic cruises in the Dixon Entrance System 208 B.3 Eulerian current measurements in the Dixon Entrance System 210 B.4 Lighthouse stations (see Fig. B.l for the geographical locations) where SSS and SST are, or have been in the past, measured with their period of observations 211 xi List of Figures 1.1 Geography of the north coast of British Columbia (Adapted from Ballantyne et al., 1996) 2.1 2 Bathymetric chart of the Dixon Entrance System. Depth contours are in metres with intervals of 50 and 500 m for regions shallower and deeper than 500 m respectively. Four regions are shown: depths less than 100 m (-), between 100 and 300 m (blank), between 300 and 500 m (+), and greater than 500 m (>) 2.2 10 Winter and summer mean sea surface atmospheric pressure maps for the North Pacific Ocean. Arrows show direction of the geostrophic winds. (Adapted from Hamilton, 1984.) 2.3 13 Mean daily discharge rate of the Skeena (dashed line) and Nass (solid line) rivers calculated for the 1928-1989 period 2.4 17 Schematic diagram of prevailing surface currents in North Pacific Ocean. (Fig. 13.17 in Thomson, 1981.) 2.5 18 Corange and cotidal values for semidiurnal tide (M ) obtained numerically. 2 The tidal range (solid line) is in cm and the tidal phase (dashed line) in degrees with Universal Time assumed. (Adapted from Foreman et al., 1993.) 21 2.6 Numerical model-derived, Eulerian tidal-residual depth-averaged current for eastern Dixon Entrance and northern Hecate Strait. (Adapted from Bowman et al., 1992.) 23 xu 2.7 Sketch of low-frequency upper water column (solid arrows) and deep (dotted arrows) currents in Dixon Entrance and northern Hecate Strait, interpreted from various measurements, model predictions and theory. (Fig. 12 in Bowman et al., 1992.) 2.8 24 Average currents and winds for the period 28 October 1984 to 1 March 1985. Solid arrows represent currents at meters within 50 m of the surface, dashes represent currents at intermediate depths, dots represent currents within 15 m of the bottom, and wide arrows represent wind at lighthouse stations. (Fig. 4 in Bowman et al., 1992.) 2.9 25 Path of a satellite-tracked drifter, drogued at 10 m, during October- December 1983. (Fig. 4 in Crawford and Greisman, 1987.) 26 2.10 Smoothed drifter tracks (tidal currents averaged out) from 1984 measurements. The number within the diamonds represent the drifter number; solid circles show position of drifters at the onset of a storm on 21 June, with winds of 15 m s _1 from the southeast. Launch and recovery dates are noted at start and end of each track. (Adapted from Crawford and Greisman, 1987.) 27 2.11 Average currents and winds for the period 22 April-19 August 1984. Solid arrows represent currents at meters within 50 m of the surface, dashes represent currents at intermediate depths, dots represent currents within 15 m of the bottom, and the wide arrow represents wind at a lighthouse station. (Fig. 5 in Crawford and Greisman, 1987.) 28 2.12 Mean daily sea-surface temperature (dashed line) and salinity (solid line), calculated for the 1939-1970 period, at Triple Island station (see Fig. 2.8). 31 2.13 Lagrangian residual current vectors, averaged over one tidal cycle, in Dixon Entrance derived from the Hecate Model (from Bell, 1963) xiii 34 4.1 3D mesh in (a) physical (x,y,z) and (b) mapped x,y,<r) coordinate space. 63 4.2 Coarse G r i d - Horizontal triangular mesh composed of 1742 elements and f 1038 nodes. Numbers represent the different open boundaries 4.3 64 Fine Grid-Horizontal triangular mesh composed of 5432 elements and 2982 nodes. Numbers represent the different open boundaries 4.4 65 Extended Grid-Horizontal triangular mesh composed of 6128 elements and 3350 nodes. Numbers represent the different open boundaries 4.5 67 Large Grid-Horizontal triangular mesh composed of 13,133 elements and 7575 nodes 4.6 68 Vertical grid with a z-level coordinate system in physical (x,y,z) coordinate space 4.7 69 Vertical linear meshes, in physical (x,y,z) coordinate space, with (a) 10 elements and 11 nodes and (b) 20 elements and 21 nodes 4.8 Depth-averaged nonlinear tidal forcing term (T (x,y,z)) rei) tum equation. Units are cm s -2 70 in the momen- and values were calculated on the Large Grid. 4.9 78 Nonlinear tidal forcing term (Q s) in the continuity equation. Units are re m s 2 -1 and values were calculated on the Large Grid. 79 4.10 Location (•) of CTD stations for cruises carried on in July 1-6 1985 (top), June 24-July 7 1991 (middle), and July 5-August 18 1991 (bottom). . . 82 4.11 Location (•) of CTD stations for cruises carried on in September 19October 10 1962 (top) and October 14-27 1984 (bottom) 83 4.12 Location (•) of CTD stations for cruises carried on in April 11-22 1984 (top) and December 2-18 1991 (bottom) 84 4.13 Horizontal homogeneous isotropic covariance model for density variations used in the statistical interpolation process xiv 87 4.14 Density field, at 10-m depth, obtained from observations and statistical interpolation for June 24-July 7 1991 (top) and December 2-18 1991 (bottom). The contours are in cr units (tr=p-1000 kg m~ ) with intervals of 3 t t 0.1 ot 88 4.15 Density field, at 40-m depth, obtained from observations and statistical interpolation for June 24-July 7 1991 (top) and December 2-18 1991 (bottom). The contours are in <r units (ir=p-1000 kg m ) with intervals of -3 t t 0.05 o- 89 t 4.16 Density field, at 100-m depth, obtained from observations and statistical interpolation for June 24-July 7 1991 (top) and December 2-18 1991 (bottom). The contours are in tr units (cr =/j-1000 kg m ) with intervals of -3 4 t 0.05 <r 90 t —* 4.17 Depth-averaged baroclinic forcing term (iZ ») for June 24-July 7 1991. Pe Units are cm s -2 and values were calculated on the Extended Grid. . . . . 91 4.18 Depth-averaged baroclinic forcing term (Rr ) for December 2-18 1991. ea Units are cm s -2 and values were calculated on the Extended Grid. . . . . 4.19 Stick vector diagrams of the hourly wind stress (H ^ ) m ret 92 for June 24-July 7 1991 at (a) buoy 145, (b) buoy 183, and (c) buoy 205 shown on Fig. B.l. Units are 10 -5 m s 2 -2 with positive values for wind blowing towards the North-West and negative towards the South-East 95 4.20 Stick vector diagrams of the hourly wind stress (H 4i ) m rea for December 2- 18 1991 at (a) buoy 145, (b) buoy 183, and (c) buoy 205 shown on Fig. B.l. Units are 10~ m s~ with positive values for wind blowing towards the 5 2 2 North-West and negative towards the South-East 4.21 Wind forcing term (v&re«) for December 2-18 1991. Units are cm s values were calculated on the Extended Grid. xv 96 -2 and 97 4.22 Daily discharge rate of the Skeena (thin line) and the Nass (thick line) rivers for the year 1991. The vertical dashed lines show June 24, July 7, December 2, and December 18 5.1 99 Residual surface elevation generated by tidal forcing. The thick solid line represents the zero surface elevation contour, thin solid lines positive values and dotted-dashed lines negative values. Contour interval is 0.2 cm. . . . 105 5.2 Residual depth-averaged flow driven by tides. Units are cm s -1 and the box delimits the zoom region shown on Fig. 5.3 5.3 Residual depth-averagedflowdriven by tides for the zoom region (box shown on Fig. 5.2). The units are cm s 5.4 106 107 _1 The five geographical regions over which the terms of the vorticity equation were evaluated, with depth contours as shown on Fig. 2.1. Regions 1 to 5 correspond respectively to Rose Spit Sill, Learmonth Bank, Langara Island, Cape Chacon, and Cape Muzon 5.5 108 Difference between the depth-averaged nonlinear tidal forcing term in the momentum equation when vertical advection is respectively neglected and —* —* —* included, T ,(x,y)-T (x,y, re (bottom). Units are cm s 5.6 z), res -2 (top) and T (x,y, ret z), zoom of Fig. 4.8, and values were calculated on the Large Grid. 113 Difference between the residual depth-averagedflowdriven by tides when neglecting the vertical advection and the residual depth-averaged flow driven by tides with vertical advection included. The units are cm s . . 115 Steady-state currents, at 10-m depth , driven by a SE wind of 10 m s . _1 5.7 _1 The units are cm s 5.8 118 _1 Steady-state currents, at 10-m depth, driven by a easterly wind of 10 m s . _1 The units are cm s . . 119 _1 xvi 5.9 Depth-averaged wind currents for a SE wind obtained by Hannah (1992). 121 5.10 Depth-averaged wind currents for a SE wind of 10 m s . The units are -1 cm s 122 _1 5.11 Current vector field, at 10-m depth, driven by the average wind observed in the December 2-18 1991 time interval. The units are cm s 123 _1 5.12 Depth-averaged tidal residual velocities on a non-rotating earth. The units are cm s 125 _1 5.13 Depth-averaged velocity field for a SE wind of 10 m s earth. The units are cm s _1 and a non-rotating 126 _1 5.14 Depth-averaged wind-drivenflowfor a SE wind and a non-rotating earth obtained by Hannah (1992) 128 5.15 Difference between the 10-m depth current vector field, driven by the average wind measured over December 2-18 1991, computed on the Grid and on the Extended Large Grid. 130 5.16 Difference between the depth-averaged tidal residual currents computed on the Fine and the Coarse grids. Units are cm s 131 _1 5.17 Tidal residual velocities in Dixon Entrance obtained by Ballantyne et al. (1996). Single shafts in multi-shafted vectors represent 5 cm s _1 132 5.18 Difference between the depth-averaged tidal residual currents computed with 11 and 21 vertical sigma-levels. Units are cm s -1 133 5.19 Difference between the depth-averaged tidal residual currents calculated with C7 o=0.005 and C o=0.01. Units are cm s 10 _1 1O 135 5.20 Difference between the depth-averaged tidal residual currents calculated with Ci—0.5 and Ci=0.3. Units are cm s xvu _1 136 6.1 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period June 24-July 7 1991. Units are cm s 6.2 -1 and the box delimits the zoom region shown on Fig. 6.2. 146 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991 for the zoom region (box shown on Fig. 6.1). Units are cm s -1 and the two circles show the approximate position of the cyclonic eddies discussed in the text 6.3 147 Current velocity field, at 40-m depth, driven by baroclinic forcing computed from the density field measured in the period June 24-July 7 1991. The units are cm s 6.4 148 _1 Current velocity field, at 100-m depth, driven by baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991. The units are cm s 6.5 149 -1 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period July 1-6 1985. Units are cm s _1 and the two circles show the approximate position of the cy- clonic eddies discussed in the text 6.6 150 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period July 5-August 18 1991. Units are cm s -1 and the two circles show the approximate position of the cyclonic eddies discussed in the text 6.7 151 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period September 19-October 10 1962. Units are cm s _1 and the three circles show the approximate position of the cyclonic eddies discussed in the text xviii 154 6.8 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period October 14-27 1984The units are cm s . 155 _1 6.9 Distribution of dynamic height (m), between 0 and 125 m, for the period September 19-October 10 1962 (Fig. 28 in Crean, 1967). Contour interval is 0.01 m (1 cm) 156 6.10 Sea-surface elevation generated by baroclinic forcing from the density field measured in the period September 19-October 10 1962. The thick solid line represents the zero surface elevation contour, thin solid lines positive values and dotted-dashed lines negative values. Contour interval is 1 cm 157 6.11 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period April 11-22 1984Units are cm s -1 and the circle shows the approximate position of the cyclonic eddy discussed in the text 158 6.12 Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period December 2-18 1991. Units are cm s _1 and the two circles show the approximate position of the cyclonic eddies discussed in the text 160 6.13 Current velocity field, at 10-m depth, driven by baroclinic forcing on a non-rotating earth. Units are cm s _1 and the forcing term was calculated from the densityfieldobserved in the period June 24-July 7 1991 xix 162 7.1 Current velocity field, at 10-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991. Units are cm s -1 and the two circles show the approximate position of the Learmonth Bank Gyre (left) and the Rose Spit Eddy (right) 7.2 166 Current velocity field, at 10-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period October 14-271984- Units are cm s _1 and the circle shows the approximate position of the Learmonth Bank Gyre 7.3 167 Locations (A,B,C) where current measurements were taken in the period September-October 1962 7.4 169 Average current components at the three locations shown on Fig. 7.3 (Fig. 27 in Crean, 1967) 7.5 169 Computed velocity profiles at the three stations shown on Fig. 7.3. The model was driven by tides, river runoff, and baroclinic forcing arising from the density field measured in the period September 19-October 10 1962. The dashed and solid lines respectively represent the west-east and southnorth components 7.6 170 Current velocity field, at 50-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period April 11-22 1984- The observed average currents are represented by bold arrows at the center of either a single, double, or triple circle and units are cm s . Note that current observations shown in double and triple circles -1 are referred to in the text. 173 xx 7.7 Current velocity field, at 150-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period 11-22 1984- April The observed average currents are represented by bold * arrows at the center of either a single, double, or triple circle and units are cm s . Note that current observations shown in double and triple circles _1 are referred to in the text 7.8 174 Current velocity field, at 50-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period October 14-27 1984- The observed average currents are represented by bold arrows at the center of either a single or double circle and units are cm s . Note that current observations shown in a double circle are _1 referred to in the text 7.9 175 Velocities, at 50-m depth, computed by Ballantyne et al. (1996). Their model was forced by tides and baroclinic pressure gradients calculated from the October 14-27 1984 cruise. Observed currents are shown in bolder type with a small circle at their foot and the mooring name nearby. Single shafts in multi-shafted vectors represent 10 cm s . (Fig. 9(f) in Ballantyne et -1 al., 1996) 177 7.10 Average currents as calculated from the summer 1991 drifter positions at 3-m depth. Estimated standard deviations are shown as ellipses around each vector. Lack of an ellipse indicates that only one observation was available (Fig. 10(a) in Ballantyne et al., 1996) 178 7.11 Current velocity field, at 3-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period July 5-August 18 1991. The units are cm s xxi _1 179 7.12 Velocities, at 3-m depth, computed by Ballantyne et al. (1996). Their model was forced by tides and baroclinic pressure gradients calculated from the July 5-August 18 1991 cruise. Single shafts in multi-shafted vectors represent 10 cm s . (Fig. 10(b) in Ballantyne et al., 1996) . . . . _1 180 B.l Locations of land wind station (X), ocean buoys (filled triangle), and lighthouse stations (•) 204 E.l Example (a) of the vertical grid and (b) of an horizontal triangular element.233 xxii Acknowledgements Je dedie cette these a la memoire de mes grands-parents Jeanne, Rene, Lydia, et Lorenzo. I would first like to thank my supervisor, Dr. Paul LeBlond, for his remarkable generosity and patience, as well as his financial support. I am very grateful to him for having given me the opportunity to attend many interesting conferences and one summer school. A special thank to Dr. Mike Foreman, who answered my numerous questions about numerical modelling and allowed me to use many of his computer codes, and to Dr. Steve Pond, who helped me understanding the physical aspects of my work. I also wish to express my gratitude to Drs Bill Crawford and Howard Freeland, who helped me in different ways, and to Roy Walters for the use of his model code. I gratefully acknowledge the financial support, through graduate fellowships, of Le Fonds FCAR from the government of Quebec as well as of the National Science and Engineering Research Council of Canada (NSERC). Un grand merci a Denis Laplante pour m'avoir maintes fois aide a solutionner mes problemes informatiques. I would also like to thank Anne Ballantyne for her help and Charles Hannah for his helpful suggestions. Thanks to Mike Woodward, who taught me about oceanographic measurements during the 3-week cruise in northern British Columbia waters, and to Bodo, Dennis, Neil, Allan, Darren, Larry, and David for their patience and the great time spent on the Tully. Thank you to Lisa, Ana and all the others I am not naming here, by fear to forget someone, for their precious friendship. Merci a mes parents, Huguette et Real, pour m'avoir toujours soutenu dans mes projets, peu importe les sacrifices que cela leur exigeaient. Grace a eux, j'ai pu realiser plusieurs de mes reves, ce doctorat en etant un. Un petit clin d'oeil a mon frere Daniel qui a toujours eu une pensee pour moi. Finalement, sans le support de mon compagnon de vie, Daniel, je n'aurais probablement jamais termine ce doctorat. Je le remercie d'avoir toujours ete la lorsque j'avais besoin de reconfort et d'encouragement. Sa grande comprehension m'a ete indispensable. Je ne pourrais terminer ces remerciements sans mentionner le petit etre qui s'est developpe en moi durant la periode de redaction de cette these. L'arrivee de mon premier enfant, Carmen, m'a procure la source d'energie necessaire pour le sprint final. XXlll Chapter 1 Introduction As field work is becoming more expensive, scientific budgets are tightening, and computer power is getting more accessible, interest in numerical modelling is increasing. Numerical modelling provides us with an inexpensive way of understanding the physical processes specific to an oceanographic region and of predicting the water currents in remote areas. However, the necessity for observations remains. In order to model the water circulation of a domain, information on the physical system is required. Observations are needed to drive the model as well as to evaluate the model accuracy. This thesis treats modelling of the mean coastal water circulation. The geographical region of interest is located on the west coast of Canada, i.e. along the north coast of British Columbia (BC). Specifically, the study area lies between 53° and 55.5°N latitude and between 129.5° and 134°W longitude. Dixon Entrance (DE) with Hecate Strait (HS), Clarence Strait, and Chatham Sound (CS) constitute the marine region studied here (Fig. 1.1). A harmonicfiniteelement model is used. The main advantage of afiniteelement model over a finite difference model is the ease with which nonrectangular geometries can be accommodated. The study region is characterized by a complex topography and bathymetry as well as rapid water property variations. A diagnostic picture of the mean circulation, responding to tidal, atmospheric, and buoyancy forcings, is studied. This thesis constitutes one step in the modelling process of these coastal waters. The development of a prognosticfiniteelement model for the region would be a follow-up to 1 Chapter 1. 136°W 2 Introduction 134° 132° 130° 128° 126° Figure 1.1: Geography of the north coast of British Columbia (Adapted from Ballantyne et al., 1996). Chapter 1. Introduction 3 the thesis. 1.1 Motivation 1.1.1 Fisheries Due to the relatively narrow continental shelf on the BC coast, most fishing operates within a zone 100 km from shore. The coastal waters along the northern BC coast are an important area for the commercial, recreational, and Native food fisheries (Petro-Canada, 1983). Crabs, coho and chinook salmon, Pacific halibut, as well as Dover sole and herring are among the commercially important fish and shellfish of the region. Large amounts of pink and sockeye salmons are also taken. The major spawning streams for salmon are the Nass and Skeena rivers (Fig. 1.1). The Skeena River ranks second as a producer of salmon in British Columbia (Hoos, 1975). Water temperature and salinity, as well as water currents and circulation, are among the important biophysical factors for fisheries and aquaculture (Ricker and McDonald, 1992). For the North Coast, water temperature is believed to be the most important factor. However, salinityfluctuationscould be a severe limitation for aquaculture, especially in Chatham Sound where there is much dilution from spring runoff discharges of the Nass and Skeena rivers. Links between physical oceanography and fisheries of this region have been established. In DE, two important water circulation features are believed to be related to fisheries: the freshwaterflushing(Crean, 1967) and the Rose Spit Eddy (Jamieson, 1985). In HS, the strong wind-driven currents are believed to carry eggs and larvae out of the strait (Ketchen, 1956; Tyler, 1986; Crawford et al., 1990). The growing awareness of the influence of oceanographic conditions onfisheriesallowed the creation, in the early 1990's, of the Ocean Production Enhancement Network (OPEN) program. This thesis work started under the OPEN program. Chapter 1.1.2 1. Introduction 4 O i l exploration Since 1972, there has been a moratorium on exploration for offshore oil and gas on Canada's West Coast. This moratorium has prevented an evaluation of this region's hydrocarbon potential (Petro-Canada, 1983). From 1967 to 1969, Shell Canada Resources Limited did exploratory offshore drilling off BC. They did not find economic amounts of hydrocarbons. Shortly after, in 1971, petroleum exploratory drilling was suspended indefinitely in response to public concern regarding the oil spill potential of tanker traffic. New geological theories and the re-evaluation of Shell's well data yield more promising interpretations of the area's hydrocarbon potential. The National Energy, Mines and Resources now encourages the exploration of Canada's domestic reserves of oil and natural gas, including those of the BC offshore. Chevron Canada Resources Ltd. (1982) acquired Shell Canada's oil and gas exploration permits in the Hecate Strait area and is seeking to renew exploratory activity in HS and Queen Charlotte Sound (Fig.1.1). Petro-Canada (1983) acquired permits for northern HS and DE. A number of hazards (e.g. accidental oil/gas blowout) that can be encountered during offshore drilling operations could result in the release of hydrocarbons into the environment. In some circumstances, the best response to an oil spill would be to monitor the movement of the slick. As oil slicks move and spread on the sea surface, the factors governing its path, or trajectory, are surface currents, winds and the Earth's rotation. To control oil spills adequately, through slick trajectory modelling, more information is needed on the strength, direction and variability of surface currents. The concern about oil spills has led to studies financed by the Panel of Energy Research and Development (PERD) on the surface circulation. A barotropic wind-driven model (Hannah et al., 1991; Hannah, 1992), a barotropic tidal model (Foreman et al., 1993), a barotropic tidal model with bottom boundary layer (Ballantyne et al., 1996), Chapter 1. Introduction 5 as well as a baroclinic diagnostic model (Ballantyne et al., 1996) have been developed for the region. This thesis constitutes a follow-up to the modelling effort started under PERD. 1.1.3 Pollution discharges Currently, the city of Prince Rupert (Fig. 1.1) discharges municipal sewage, street sewer, storm runoff water, and commercial-industrial effluents into the harbour (Stucchi and Orr, 1993). This practice of sewage disposal to the harbour has created high conform counts near shore and is aesthetically unpleasant. The city is now proposing a sewerage extension program. On top of that, Prince Rupert is the proposed support base for Chevron's activities and would be the most likely location for the shore base of PetroCanada's offshore operations. In the Prince Rupert Harbour complex, a study is required in order to define its complex circulation patterns through all stages of tidal change and to predict the fate of its ever-increasing pollution discharges. Because the harbour opens directly onto Chatham Sound, the oceanography of the sound directly influences that of the harbour. CS provides the source water for Prince Rupert Harbour. A circulation model could be useful to track the dispersal and movement of effluent and contaminants introduced to the system as well as to optimize locations and configurations of effluent discharge. A circulation model for a greater region would help in the specification of boundary conditions in local models of effluent diffusion. 1.2 Objectives One of the main thesis objectives is to increase our understanding of the dynamics that control the mean water circulation along the northern BC coast. This region is known to be highly baroclinic and thus barotropic models can not adequately represent the Chapter 1. Introduction 6 physical system. Emphasis is put here on modelling the tidal residuals with stratification effects, the tri-dimensional (3D) wind-driven circulation, and the buoyancy-driven coastal circulation. The baroclinic model results are compared with results of earlier models and existing observations. Another goal is to characterize the regional mean water circulation for different times of the year. The improvement of coastal modelling for complex regions is the third important purpose. To achieve these goals, numerous modelling steps are necessarily: • Incorporation of the essential physics into the governing equations • Creation of a spatial grid used to discretize the governing equations • Choice of friction parameterization • Choice of appropriate open boundary conditions applied to the model domain • Preparation of the observations used to drive the model • Use of an objective analysis technique in order to interpolate density data such that baroclinic pressure gradients can be calculated • Evaluation of model results using Eulerian and Lagrangian current observations 1.3 Overview of Thesis The thesis is organized as follows: Chapter 2 contains a review of the study area's physical environment. The geography, sea-floor sedimentology, meteorology and climatology, hydrology, as well as physical oceanography and climatic change are discussed. A review of laboratory and numerical models developed in the past for this region is also provided. Chapter 1. Introduction 7 Chapter 3 introduces the physics of the model as well as the numerical solution procedure. In chapter 4, the essential steps for modelling the read physical system are presented. Spatial discretization of the domain, open boundaries, model parameters, and model forcings Eire discussed. The chapter 5 deals with tidal rectification and wind-driven circulation. Topics include the importance of vertical advection of tidal momentum, the influence of the earth's rotation, the effect of stratification, and the sensitivity of the model. Model results are compared with earlier studies. Chapter 6 gives model results obtained with river runoff and baroclinic forcings. In chapter 7, the model ability to reproduce the real physical system is evaluated. The complete model solution, currents driven by all forcings (tides, wind, river runoff, and baroclinic pressure gradients), is compared with current observations. Chapter 8 concludes the thesis with a summary of the research work and suggestions for future investigations. The first appendix is a summary of the diverse symbols used through the thesis. A historical review of the observational programmes in the region of interest is presented in the second appendix. Mathematical details of the model and of the statistical interpolation used are given in the next four appendices. Chapter 2 R e v i e w of the Physical Environment and Earlier M o d e l s This chapter constitutes an overview of the physical environment as known, to the best of my knowledge, from observations (see Appendix B) and earlier modelling efforts. Five main aspects are discussed: regional topography and bathymetry; sea floor sedimentology; meteorology and climatology; hydrology and physical oceanography. Climatic variability of this environment is also discussed. The chapter ends with an overview of the laboratory and numerical models developed to study the physical oceanography of the region studied in this thesis. The main pioneering works and results have been reported by Trites (1956) and Crean (1967). Dodimead (1980) published a review which included the work of Crean but not that of Trites. Shortly after, Thomson (1981) reviewed the physical environment of the region for the general public. A summary, which includes research results obtained up to the 1980's, was more recently presented by Thomson (1989). Good reviews of the physical environment, as well as of the biological and socioeconomic environments, of the region of interest can also be found in Chevron Canada Resources (1982), Petro-Canada (1983), and Pucker and McDonald (1992). The Skeena Estuary review by Hoos (1975) is a good document as well but deals only with a small subsection of the study area. 2.1 Topography and Bathymetry The marine region of interest, referred through the thesis as the Dixon Entrance System, is the body of water which separates Graham Island, the northern section of the Queen 8 Chapter 2. Review of the Physical Environment and Earlier Models 9 Charlotte (QC) Islands, from the BC mainland and the Alaska Panhandle (Fig. 1.1). It includes Dixon Entrance, its two main adjacent seaways, Hecate and Clarence Straits, as well as an estuary called Chatham Sound. The regional topography is characterized by a high relief along the mainland coast where rugged coastal mountains rise up to 2000 m and contain numerous fjords and inlets. On islands off the mainland coast, the relief is low, generally below 200 m in elevation. On Graham Island, the eastern side is relatively low lying and the western part mountainous. An interesting feature of the Pacific continental shelf, to the west of QC Islands, is its remarkable narrowness; the 1000 m depth contour lies within only 25 km of most of the coastline. To the north, along Graham Island, the shelf widens and reaches a maximum width of roughly 30 km at 54° N latitude. There it merges with the broad continental shelf west of Dixon Entrance and the continental slope becomes considerably less steep. Further north, the shelf widens to approximately 100 km. 2.1.1 Dixon Entrance and adjacent seaways Dixon Entrance has a width of about 55-65 km and a length of approximately 150 km. This broad east-west depression in the continental shelf is bounded to the north by two mountainous Alaskan islands, ranging from 300-1500 m in elevation, and to the south by Graham Island. DE is indeed a gap between two mountain barriers. This westward deepening depression reaches, near the seaward entrance, a depth of 400 m. The deep entrance channel is however broken by a broad shallow bank into North and South entrances (Fig. 2.1). The bank, called Learmonth Bank, is about 19 km in length, oriented northwest (NW) to southeast (SE), has a width of about 8 km, and rises to within 35 m to the surface. To the east, DE shoals gradually until it reaches a sill (270 m deep) south of Cape Chacon (Fig. 1.1). From there, depths increase as the channel swings 10 Chapter 2. Review of the Physical Environment and Earlier Models 134 W 134 W 133 W 133 W 132 W 132 W 131 W 131 W 130 W 130 W Figure 2.1: Bathymetric chart of the Dixon Entrance System. Depth contours are in metres with intervals of 50 and 500 m for regions shallower and deeper than 500 m respectively. Four regions are shown: depths less than 100 m ( - ) , between 100 and 300 m (blank), between 300 and 500 m (+), and greater than 500 m (>). Chapter 2. Review of the Physical Environment and EarHer Models 11 northward into a deep (400 m) trench named Clarence Strait. The latter, a 200 km long seaway, narrows and shoals to the north. A broad and shallow (20-40 m deep) bank, labelled Rose Spit Sill, separates the southeastern part of Dixon Entrance from the north part of Hecate Strait (Fig. 2.1). This strait is a strong asymmetrical trough, about 220 km long, opening southward. It is characterized by a gradual increase in cross-sectional area from north (about 60 km wide and 20 to 80 m deep) to south (approximately 130 km wide and 300 m deep) (Crawford et al., 1988). The main axis of HS, a narrow long SE-NW submarine valley, hugs the BC mainlandflank.On its western side, the strait shallows rapidly to form a wide and shallow (10-60 m deep) platform which separates the main axis from the QC Islands. 2.1.2 Chatham Sound Chatham Sound, a semi-enclosed basin, is separated from both DE and HS by a submarine ridge supporting a chain of shoals and islands. At its northern end, the ridge is broken by a deep, irregular channel which extends westward to join Clarence Strait (Fig. 2.1). Depths in the central and northern parts of the basin are usually in excess of 100 m. The southeast part is shallower (50 m deep) because of millennia of sediment addition by the Skeena River. Unlike most other estuaries, there is no distinct delta formed by this river. Rather, there is a series of pockets of shallow water where sediments have built up to form shoals and sand banks. Part of this pattern is due to the strong tidal currents that pass through the many passages, preventing sedimentation along the broad delta-front of the Skeena River estuary. Due to shoals and sand banks, CS is actually one of the most difficult areas to navigate along the BC coast; sizes, shapes, and locations of shoals and banks can change rather rapidly (Hoos, 1975). These banks are however important since they provide nursery areas for juvenilefishspecies. Another major river, the Nass, alsoflowsinto the sound. In contrast to the Skeena, this river does not flow Chapter 2. Review of the Physical Environment and Earlier Models 12 directly into CS but rather via an inlet (Portland Inlet). Therefore, the sediments from the Nass do not settle in the sound but in the inlet which is deep enough (300-400 m) to trap any settling material (Crawford and Thomson, 1991). Due to the introduction of fresh water from two important separate sources, to the great water depth variations, and to its very irregular coastline, Chatham Sound constitutes a complex estuary. 2.2 Sea F l o o r Sediments In Dixon Entrance, theflooris generally covered with a relatively thin veneer of pebbles to boulders, with or without sand or mud, and with occasional outcropping bedrock (Mathews, 1958). The north side and main axis of DE, even in its deepest part, are covered by clean gravels, suggesting that theflooris scoured by bottom currents. In contrast, the area just north of Graham Island, on the south side of DE, consists of sand, gravel, and boulders. A fine sand, often silt, slope extends northeastward off Rose Spit (Fig. 1.1) and becomes dominantly mud at about the 100 m-depth line. From there, it gives way to gravel and rock before reaching Dundas Island (Fig. 1.1). Erosion of unconsolidated sediments from Graham Island is the origin of the sand beach that extends up to Rose Spit and which is still being constructed progressively northeastward. On the western half of Hecate Strait, i.e. the platform off QC Islands,finematerials are mostly absent and sand and gravel predominate (Westrheim et al., 1980). On the eastern side of the strait, more specifically in lower parts of the axial trough (below 100 m), mud, often mixed with sand, is common. The surficial geology of Dixon Entrance, Hecate Strait, and the shelf area to the west of the QC Islands is still relatively poorly known. "These regions are considerably more complex than the southern shelves in terms of roles played by active tectonics, wave and Chapter 2. Review of the Physical Environment and Earlier Models 13 Figure 2.2: Winter and summer mean sea surface atmospheric pressure maps for the North Pacific Ocean. Arrows show direction of the geostrophic winds. (Adapted from Hamilton, 1984.) tidal regimes, and coastal erosion and sediment transport" (Bornhold and Barrie, 1991). 2.3 Meteorology and Climatology The climate of the Canadian west coast is usually mild, with cool summers and warm, wet winters. The prevailing onshore flow of marine air, whose temperature is regulated by heat exchange with the upper ocean, is the cause of that climate. Two semi-permanent atmospheric pressure systems characterize the North Pacific Ocean (Fig. 2.2). The Aleutian Low, centered over the Bering Sea and Gulf of Alaska, dominates in winter while the North Pacific High, centered off California, dominates in summer. During the winter, storms are frequent in northern B C coastal waters (Kendrew and Kerr, 1955). Chapter 2. Review of the Physical Environment and Earlier Models 2.3.1 14 Winds Prevailing winds in the northeast Pacific are controlled by the location and intensity of the two major atmospheric pressure systems. During the winter, the Aleutian Low produces strong south-to-southeasterly winds along the coast. In summer, north-to-northwesterly winds are produced by the North Pacific High. Westward migrating synoptic-scale atmospheric systems (highs and lows) are responsible for modification, over period of days to weeks, of these prevailing winds. Synoptic winds are modified in near-coastal areas by the mountainous terrain. Within 50 km of shore, winds are constrained to blow parallel to the coast. Therefore, prevailing wind directions tend to be aligned west-east and NW-SE in Dixon Entrance and Hecate Strait respectively. During winter (October-March), dominant winds are from the east in DE and the southeast in HS. Note that the strongest winds are then found in central HS. In DE, winter winds are usually slightly stronger and persist longer at the eastern end, which is exposed to SE gales from HS, than at the western end. In summer (AprilSeptember) , winds tend to be weaker and respectively more westerly and northwestherly in DE and HS. The strongest winds occur in November through February, and weakest winds in July and August. Late September to early October usually marks the time of an abrupt increase in mean wind speeds. Significant local variations can arise from other effects. For example, north and northeasterly outflow of Squamish-type winds, essentially an atmospheric density-current, often occur within coastal areas in winter. Also, the near-coast winds can be modified by the diurnal sea-breeze set up by daily heating and cooling of the adjacent land mass. Chapter 2. Review of the Physical Environment and Earlier Models 2.3.2 15 Precipitation Heavy precipitation, due mainly to the upward motion of air caused by frequent travelling storms, characterizes the entire northern British Columbia coast. The annual precipitation, spatially averaged, reaches 254 cm (Royer, 1982). The range goes from less than 160 cm over eastern QC Islands to more than 320 cm along the western slopes of the coastal mountains. In general, the amount of precipitation increases from south to north and from the interior to the coastal ranges. Through the course of a year, there is a marked variation in the precipitation rate. In winter, the Aleutian Low creates storms with high moisture content. The coastal mountain ranges, bordering the Gulf of Alaska, act as a barrier to these storms which move eastwards across the North Pacific. Adiabatic elevation of these moist air masses over BC mountains cause very high rates of precipitation in the form of rain and snow. In summer, the high pressure system obstructs the eastward flow of moisture-laden air and produces a summer precipitation minimum. With the onset of fall, precipitation increases rapidly such that a maximum is reached during the October-December period. As from November to March the bulk of the precipitation is generally retained as snow (Kendrew and Kerr, 1955), almost down to sea level, the maximum direct freshwater contribution is released into the coastal seaways in October. 2.4 Hydrology The water released by the atmosphere over the coastal region is either lost by evaporation, discharged as runoff within a year, or held for longer periods by local glaciers. Evaporation is a minor process for British Columbia (Royer, 1979). As regards the local glaciers, they provide storage and release of water with a residence time that can reach up to thousands of years (Fountain and Tangborn, 1985). As most of the basins in northern Chapter 2. Review of the Physical Environment and Earlier Models 16 BC are glacierized, glaciers have a local effect on runoff. The annual hydrologic discharge pattern for streams on the North BC coast is a complex one with various factors altering the timing of maximum runoff. Coastal streams on low relief regions usually have peak runoff in autumn to early winter. This is a consequence of the arrival of the Aleutian Low pressure cell and associated heavy precipitation events in the form of rain. Coastal streams in area of high relief have two peak runoff periods during the year, one in autumn and one in spring. The autumn peak is in response to the Aleutian Low effect while the spring one is due to the melting of the stored snowpack, which has built up over the course of winter and early spring. The relative magnitude of the spring and autumn maxima depends not only on the location and average elevation of the drainage basin, but also on the annual cycle of precipitation appropriate to the particular area. 2.4.1 Skeena and Nass rivers Two main rivers, the Skeena and the Nass, flow into northern BC waters. The former empties directly into Chatham Sound, 20 km south of Prince Rupert. The latter discharges into the sound but through Portland Inlet which has, through its lengthy tributary fjords, considerable direct glacier run-off into it. The drainage basin dimensions are, for the Skeena and Nass rivers, 42 200 km and 19 200 km respectively. The basins 2 2 are topographically dominated by mountains. As the time of maximum high elevation snow melt from those mountains is usually the month of June, the peak flow of both rivers is in June. The rivers' hydrographs are thus dominated by a massive spring runoff event, beginning in late April with maxima in the first weeks of June, followed by a gradual reduction in outflow until October, when a secondary peak associated with heavy precipitation occurs (Fig. 2.3). Note that location of the gauges some distance upstream Chapter 2. Review of the Physical Environment and Earlier Models Jan. Mar. May July Sept. 17 Nov. Figure 2.3: Mean daily discharge rate of the Skeena (dashed line) and Nass (solid line) rivers calculated for the 1928-1989 period. from the river mouth, 145 km and 75 km for the Skeena and Nass respectively, must contribute to some extent to the apparent weakness of the secondary maximum; the gauged river runoff is not being affected by much of the heavy precipitation near the coast. Note also that because autumn storms are coastal in nature, the hydrographs for the Skeena may miss very important downstream discharge events. Finally, it should be pointed out that many small local rivers flow in to the coastal waters and, in the aggregate, could also account for an important freshwater input to the northern coastal seaways (LeBlond et al., 1983). 2.5 Physical Oceanography Oceanographic conditions, i.e. water property structures and currents, along the northern B C coast are determined in part by the adjacent oceanic environment and in part by coastal factors. Among the latter are tides, winds, river discharge, coastline geometry, and bottom topography. Currents can be separated into two basic categories: tidal currents associated directly Chapter 2. Review of the Physical Environment and Earlier Models 160° w 180 F ? WESTERN SUBARCTIC GYRE — 18 BERING SEA GYRE iV PAPA Subarctic Curnnt West Wind Drift SUBARCTIC BOUNDARY North PK/flc Currant V V Figure 2.4: Schematic diagram of prevaihng surface currents in North Pacific Ocean. (Fig. 13.17 in Thomson, 1981.) with the luni-solar tide generating force, and non-tidal currents associated with all other physical processes. Included in the broad second category are the directly wind-induced coastal currents, buoyancy-driven coastal currents, and estuarine circulation generated by land-derived runoff. Tidal currents are a superposition of flows with specific periodicities (semidiurnal, diurnal, fortnightly, etc.) and can be readily extracted from the data provided the record is sufficiently long. The non-tidal circulation is what is left after the tidal contribution is removed and is not as readily determined; it is much more difficult to analyze and interpret. 2.5.1 Adjacent oceanic currents Near the west coast of North America, the Subarctic Current, an eastward flow, divides into two branches (Fig. 2.4). A northern branch curves to the northeast into the Gulf Chapter 2. Review of the Physical Environment and Earlier Models 19 of Alaska, as the Alaska Current, and a southern branch turns to the southeast as the California Current. The bifurcation region is one of confused currents within the region between 45°-50° N and 130-150° W. Off the continental shelf of the QC Islands, the Alaska Current, an eastern boundary current,flowsnorthward as a broad, variable, and slow drift. Speeds range from an average of about 25 cm s cm s _1 in winter, but may be accelerated to over 50 cm s _1 in summer to around 35 _1 by strong southeast winds in winter (Thomson, 1981). Closer to the shores of QC Islands, an intensified flow, the Haida Current, has been identified by Thomson and Emery (1986). It is a narrow (20 to 30 km wide and over 200 km long), warm currentflowingpoleward, from the northern part of Moresby Island (Fig. 1.1) to the southern portion of the Alaska Panhandle, over the continental slope. This feature appears to be particularly well established in winter and early spring. The current maintains its integrity as itflowspast the mouth of Dixon Entrance and does not, as one might expect, attempt to hug the coastline. The presence of a west-east gradient in sea surface elevation (positive eastward) could set-up a barotropic pressure gradient to counter any tendency of the Haida Current to enter the entrance (Thomson, 1989). On another hand, theflowseems to be susceptible to baroclinic instabilities once it reaches the vicinity of DE. 2.5.2 Tides and tidal currents Tides, within the region of interest, are mixed and predominantly semidiurnal. They are dominated by the principal lunar semidiurnal constituent (M ), and the principal mixed 2 diurnal constituent (Ki), and co-oscillate with tides in the adjoining North Pacific Ocean. The tide propagates northward along the Pacific coast and as it reaches the entrance of Queen Charlotte Sound, one branch enters Hecate Strait and another propagates north along the west coast of QC Islands. The former then propagates north along the mainland Chapter 2. Review of the Physical Environment and Earlier Models 20 coast while the latter travels east through Dixon Entrance. Those two branches meet near the north end of HS. The combined tides then swing westward across HS to reach the southeast shores of Graham Island. Fig. 2.5 shows the co-amplitudes and co-phases of the M tide. 2 Shoaling effects cause the mean range of the tide to increase from 3.0 m across the mouth of QC Sound to 4.8 m midway along HS, and from about 3.5 m at the mouth of DE to 5.0 m at Prince Rupert (Thomson, 1981). In Chatham Sound, the tidal range is relatively large, having a mean value of about 6.0 m (Trites, 1956). Actually, the tidal range in the vicinity of the Skeena River is the largest on the Canadian Pacific coast, exceeded only in some inlets where funnelling effect amplifies the tidal range (Hoos, 1975). Large volumes of water have thus to be moved in and out of CS, through constricted passages, and strong currents (1 to 2 knots) as well as intensive tidal turbulence are created. In Hecate Strait, tidal streams tend to move parallel to the main axis of the channel. They are basically rectilinear, due to restrictions on cross-strait flow by the valley-like bathymetry. Under such conditions, the stream simply decelerates and accelerates in the opposite direction. Floods are therefore directed to the north and ebbs to the south at maximum rates of roughly 50 cm s . _1 To the north of the meeting area of the tides, there is a reversal in the ebb andflooddirections; principal ebbs are directed to the north and principalfloodsto the southeast. Along the eastern shores of the relatively narrow northern part of HS, the ebb is stronger than the flood, suggesting a net tidal flow northward through the strait and into DE. The confluences of tides in combination with the right-angle geometry of the coastline further complicates the tidal stream within the central portion of Dixon Entrance. Instead of having near uniform strength, flood streams are decidedly stronger on the southern side of DE and ebb streams stronger on the northern side. South of Cape Chacon, where the ebb from Clarence Strait converges Chapter 2. Review of the Physical Environment and Earlier Models 21 Figure 2.5: Corange and cotidal values for semidiurnal tide (M ) obtained numerically. The tidal range (solid line) is in cm and the tidal phase (dashed line) in degrees with Universal Time assumed. (Adapted from Foreman et al., 1993.) 2 Chapter 2. Review of the Physical Environment and EarUer Models 22 with the ebb that originates from northern HS, an intense westward-setting stream is formed. Further toward the mouth of DE, tidal streams run parallel to the trend of the channel axis and eventually merge with the northward setting flood and southward setting ebb along the outer coast. Tidal currents can be made more complicated as a result of vertical density gradients, nonlinear effects, coastal geometry, or interaction of theflowwith bottom topography. Such factors can be responsible for processes as diverse as internal tides, continental shelf waves, and tidal rectification. Baroclinic tides and continental shelf waves of tidal period are presumably generated along the west coast of Graham Island through interaction of the astronomical tidal streams with the coastal bottom topography (Thomson, 1989). Tidal rectification is believed to be important along the Rose Spit Sill, south of Cape Chacon, and around Learmonth Bank (Figs 2.6 and 2.7). 2.5.3 N o n - t i d a l currents The wind and water density driven coastal currents, which vary in both time and space, are difficult to predict in the short term. The seasonal patterns only can be estimated with acceptable accuracy. In winter, the wind-driven circulation predominates. Surface water movement in Hecate Strait is a general northwardflowcaused by the predominant southeasterly winds. If the winds are strong, a current northward through HS and Clarence Strait results, with seaward (westerly) current through Dixon Entrance and then northward again along the Alaskan coast. If the southeasterly winds are weak, the current will be largely confined to Hecate and Clarence Straits. Examples of observed winter currents are shown on Figs 2.8 and 2.9. In late spring and early summer, the southeasterly winds subside and large volumes Chapter 2. Review of the Physical Environment and Earlier Models Cape Chacon scale (cm s") 1 6 0 km Figure 2.6: Numerical model-derived, Eulerian tidal-residual depth-averaged current for eastern Dixon Entrance and northern Hecate Strait. (Adapted from Bowman et al., 1992.) 23 Chapter 2. Review of the Physicsd Environment and Earlier Models 24 Figure 2.7: Sketch of low-frequency upper water column (solid arrows) and deep (dotted arrows) currents in Dixon Entrance and northern Hecate Strait, interpreted from various measurements, model predictions and theory. (Fig. 12 in Bowman et al., 1992.) Chapter 2. Review of the Physical Environment and Earlier Models 25 Figure 2.8: Average currents and winds for the period 28 October 1984 to 1 March 1985. Solid arrows represent currents at meters within 50 m of the surface, dashes represent currents at intermediate depths, dots represent currents within 15 m of the bottom, and wide arrows represent wind at lighthouse stations. (Fig. 4 in Bowman et al., 1992.) Chapter 2. Review of the Physical Environment and Earlier Models 26 135 Figure 2.9: Path of a satellite-tracked drifter, drogued at 10 m, during October- December 1983. (Fig. 4 in Crawford and Greisman, 1987.) of snow melt are discharged from the multitude of coastal estuaries. Freshwater discharges from the Nass and Skeena rivers reach their maxima (Fig. 2.3). The general northward moving currents weaken and a wide-spread estuarine circulation, northward in HS and seaward through DE, is formed. Note that runoff has only a small effect on the circulation of HS. Due to the hydraulic pressure head set up by the accumulation of fresh water along the mainland shores, an upper zone brackish water which is forced to move seaward results. As it progresses, the upper zone becomes increasingly more saline through entrainment of oceanic water from below. Observations suggest that the seaward movement of brackish water is generally confined to the northern side of DE except in the vicinity of Learmonth Bank, where topographically induced recirculation occurs (Thomson and Emery, 1986). A strong outflow past Langara Island, which is on the left side of the entrance, looking seaward, has been observed with drifters (Fig. 2.10). Below the surface layer, a reversed landward (eastward) intrusion of cool saline water, considerably slower than the surface outflow, is created to compensate the outflow. The 27 Chapter 2. Review of the Physical Environment and Earher Models Figure 2.10: Smoothed drifter tracks (tidal currents averaged out) from 1984 measurements. The number within the diamonds represent the drifter number; solid circles show position of drifters at the onset of a storm on 21 June, with winds of 15 m s from the southeast. Launch and recovery dates are noted at start and end of each track. (Adapted from Crawford and Greisman, 1987.) - 1 Chapter 2. Review of the Physical Environment and Earher Models 28 Figure 2.11: Average currents and winds for the period 22 April-19 August 1984. Solid arrows represent currents at meters within 50 m of the surface, dashes represent currents at intermediate depths, dots represent currents within 15 m of the bottom, and the wide arrow represents wind at a lighthouse station. (Fig. 5 in Crawford and Greisman, 1987.) inflow exists at depths deeper than approximately 50 m except along the northern side of Learmonth Bank, where a landwardflowexists at all depths. Also, westerly winds can develop in DE and cause a net eastward flow of water at all depths everywhere in the entrance. Examples of observed summer currents are shown on Figs 2.10 and 2.11. By late autumn, the nontidal component of the circulation becomes dominated by winds, and what remains of the brackish upper zone isflushedseaward through Clarence Strait. This condition prevails throughout the winter and early spring until the return of warmer conditions and melting of snow in the high elevations of the coastal mountains. The presence of a net northwardflowthrough HS was first observed with drift experiments made by Thompson and Van Cleve (1936). From these experiments, the seasonal movement of oceanic water into DE was also observed. Through Clarence Strait, Haight (1926) has inferred a net northwardflow.As the strait narrows and shoals to the north, any net transport through this passage must however be small (Crawford and Greisman, Chapter 2. Review of the Physical Environment and Earlier Models 29 1987). 2.5.4 Gyral circulation The presence of a cyclonic vortex, whichfillsthe entire eastern end of Dixon Entrance, wasfirstsuggested by Crean (1967) from examination of dynamic heightfieldsobtained using temperature and salinity profiles taken in 1962. Subsequentially, it became referred to as the Rose Spit Eddy (Fig. 2.7). It is only in 1984 that direct current observations of the gyre were obtained. They included results from long-term current meter moorings as well as Lagrangian drifters deployments (Crawford and Greisman, 1987) . The period of rotation then obtained for the Rose Spit Eddy was about 20-30 days. More recently, the period of rotation was estimated from drifter measurements made in 1991 to be 5-40 days, depending on distance from center of the eddy (Crawford, 1995). This eddying motion has also been observed in satellite thermal images, during periods when climatic conditions are conducive to producing thermally distinct water masses within northern Hecate Strait (Visser et al., 1990). Occasionally, the Rose Spit Eddy appears to break down, when strong winds, during a fall weather transition in October, hit the region (Bowman et al., 1992) . Even though the eddy occasionally disappeared, it is normally present and dominates the averageflow.Visser et al. (1990) showed, using a tidal stress approach, that the generation of a basin-wide barotropic eddy in eastern DE is consistent with tidal rectification over the Rose Spit Sill. A clockwise circulation around the easternflankof Learmonth Bank, in western DE, has also been observed (Crawford and Greisman, 1987). The eddy is labelled Bank Gyre Learmonth and is claimed to be driven by tidal rectification. The gyre is apparently not always present but when it is, it appears to combine with the western side of the Rose Spit Eddy to carry surface water southward across the entrance. Then, the current bifurcates; a significant component exits along the south coast of DE, into the Pacific Ocean, and Chapter 2. Review of the Physical Environment and Earher Models 30 the remainder recirculates as the southern arc of the Rose Spit Eddy. Observations do not however reveal the degree of closure of theflowwest of Learmonth Bank; the eddy may or may not be closed. Observations indicate that those gyres are confined to the upper 50-100 m of the water column (Crawford and Greisman, 1987). The eddies seem to be mainly driven by tidal rectification over steep topographic features while variations in theflowseem to be driven primarily by variations in the winds. 2.5.5 Water properties The distribution of temperature within the sea is primarily determined by heat transfer at the air-sea interface and by advective and turbulent processes beneath the sea surface. In the vicinity of Dixon Entrance, the air-sea heat exchange is an oceanic gain during the summer, due to solar radiation, and an oceanic loss throughout the remainder of the year, with the maximum occurring in January, due to other processes such as evaporation (Tabata, 1957). The sea-surface temperature (SST) is mainly dominated by the annual heating cycle (Fig. 2.12). However, relatively cold water in summer and warm water in winter, brought by wind-driven currents, partially counteract the effect of solar heating. The reverse seasonal variation in deep waters also contributes to decreasing the effect of solar heating. Coastal deep waters are colder in summer than in winter due to the summer northwest winds and subsequent upwelling. The SST can also be influenced by the seasonal effects from runoff temperatures. Water brought by the rivers acquires some of the temperature characteristics (warm in summer and cold in winter) of the province's interior during its long journey to the sea. Nevertheless, tidal mixing tends to distribute the effects over the water column and this temperature impact is rather small. Coastal sea-surface salinity (SSS) varies mainly seasonally in accordance with the input of freshwater, from both precipitation and runoff. In the region of interest, the low Chapter 2. Review of the Physical Environment and Earlier Models Jan. Mar. May July Sept. 31 Nov. Figure 2.12: Mean daily sea-surface temperature (dashed line) and salinity (solid line), calculated for the 1939-1970 period, at Triple Island station (see Fig. 2.8). salinity of inshore waters is primarily due to the Skeena and Nass rivers discharge (see Figs 2.3 and 2.12). The influence of freshwater discharge is more strongly felt in Dixon Entrance than in Hecate Strait, and in DE is more strongly felt along the northern side than the southern side. In opposition to the surface waters, deep waters are saltier in summer than in winter. The intensified estuarine circulation and northwesterly winds and subsequent upwelling during summer time are the cause. Throughout the year, water density off the northern BC coast is determined primarily by salinity rather than temperature. In winter, the vertical gradient of density is almost solely determined by salinity while in summer the main effect of surface warming is to further enhance the strong density gradient associated with the estuarine discharge. The upper mixed layer of the study region is deeper in winter than in summer. Mechanical mixing by winds with reduced stratification near the surface is the main cause. In the summer, because surface heating and freshwater runoff stratify the water column, the upper layer is thinnest. The higher SST regions are associated in the summer with areas of low-salinity. It has been suggested that this association is due to the high stability of Chapter 2. Review of the Physical Environment and EarHer Models 32 the surface waters, and consequent reduction in the downward transfer of heat by eddy diffusion (Pickard and McLeod, 1953). It is of interest to note that in shallow waters, characterized by intense tidal currents, strong vertical mixing occurs such that surface waters are mixed with colder deep water and little or no temperature stratification is present. The presence of a persistent halocline (a vertical salinity gradient) throughout the greater part of DE has been identified. A sharp halocline exists however only during the spring and summer runoff period. The annual variations of the halocline depth, minimum in summer and maximum in winter, are more pronounced in the western part of DE. The growth of the thermocline from spring to summer and its decay from autumn to winter are however similar throughout the entrance. Water column stratification in DE results in an internal Rossby radius of approximately 20 to 30 km, representing about one third (in summer) to one half (in winter) of the width of the entrance (Bowman et al., 1992). 2.6 Climatic Variability In British Columbia, an increase of 5.9 cm in the annual precipitation has been observed from the 1901-1930 to the 1931-1960 periods (Powell, 1965). The increase was found to be slightly higher along the coast than in the interior and higher in winter than in other seasons. Also, the average discharge of most BC rivers was found to have increased by about 20% in the mid-1940's and to have remained at this level for about 30 years (Slaymaker, 1990). An increase in precipitation, mainly during the winter months, was held to be responsible. The only exceptions were heavily glaciarized basins and basins in the semi-arid southern Interior Plateau. Note that the precipitationfluctuationsare a more important factor in runoff changes of the region of interest than global temperature fluctuations (Karl and Riebsame, 1989). It was actually established that runoff changes Chapter 2. Review of the Physical Environment and Earlier Models 33 in BC rivers, over the 1913-77 period, were broadly controlled by precipitation changes and, in the majority of cases, the effect was evident in the same year as the precipitation change. A growth in the extent of glaciers, starting in the 1960's and reaching a maximum extension or volume in the 1970's and early 1980's, has also been observed (Ricker and McDonald, 1992). The relatively cold and moist winters of the 1940's and 1950's are likely the cause of this expansion. Changes in surface water properties have been observed. The SST time series at Langara Island (see Fig. B.l) shows a decreasing trend of 0.53°C per century over the 1936-89 period (Freeland, 1990). At the same location, a decreasing SSS trend of 0.4 %o has also been observed over the 1936-1970 period (Webster and Farmer, 1976). 2.7 Physical Oceanographic Models 2.7.1 Laboratory model Thefirstphysical model to be developed to study the water circulation along the northern BC coast was a non-rotating hydraulic tidal model called the Hecate Model (Bell and Boston, 1962; 1963). The tides were produced by electronically-programmed tilting plates located across the mouth of Queen Charlotte Sound and of Dixon Entrance. Except in constricted coastal passages, the model predicted instantaneous tidal velocities of about 50 cm s -1 (Bell, 1963a; 1963b). Since the model was not rotating, no Coriolis force was present and the water motion could not correctly reproduce tidal current ellipses. The main feature obtained with the model is a strong net cyclonic vortex in DE, centered approximately between Cape Chacon and Rose Spit (Fig. 2.13). As winds and thermohaline circulation were absent in the experiment, the process of tidal rectification was the only mechanism left to generate the eddy. When the northern end of Hecate Strait was blocked by sandbags, that vortex was totally ehminated. It thus appeared that Chapter 2. Review of the Physiccd Environment and EarUer Models 34 Figure 2.13: Lagrangian residual current vectors, averaged over one tidal cycle, in Dixon Entrance derived from the Hecate Model (from Bell, 1963). the major factor in producing the basin-wide residual cyclonic vortex was the strong tidal movement adjacent to Rose Spit, suggesting that the strong net cyclonic vortex resulted from the meeting of the Q C Sound and D E tides in northern HS. Note that the mean currents obtained north of Rose Spit reached as high as 50 cm s , considerably in excess _1 of measured currents, typically about 6-22 cm s . Bowman et al. (1992) brought up the -1 possibility that the eddy was an artifact that could have been generated to a significant degree by phase errors. 2.7.2 Barotropic numerical models To the best of my knowledge, the first numerical model developed for the study area was by Kinney et al. (1976). The purpose was to predict possible oil spill movements for the Kitimat Pipeline Project. It was a M barotropic model, with 9.3 km resolution. 2 Subsequently, a 18-km resolution barotropic tidal model, which included both the M and K i tidal constituents, was developed by Flather (1987). 2 Chapter 2. Review of the Physical Environment and Earlier Models 35 Bowman et al. (1992) developed a barotropic tidal model in order to investigate the importance of the earth's rotation on tidally-rectified currents generated in the vicinity of the Rose Spit Sill. Simplified geometry was used for both Dixon Entrance and Hecate Strait; they were defined as two uniform channels. The most interesting feature obtained was a pronounced along-isobath current over the Rose Spit Sill (Fig. 2.6). The Coriolis effect was found to be important in producing the residual flow along that sill. A headland jet, with velocity of about 5 cm s , was also produced south of Cape Chacon. It was _1 found that the jet arose primarily from inertial effects, rather than Coriolis effects. The model showed some evidence of the Rose Spit Eddy but the period of rotation found (100150 days) was much larger than what the observations (20-30 days and 5-40 days) and the hydraulic model (10-15 days) suggest. Thus, barotropic tidal rectification alone could not explain the presence of this eddy. However, it was suggested that tidal rectification over the Rose Spit Sill and near Cape Chacon was important in recirculating some of the flow southward across, and eastward around the basin on the south side, helping to close the gyre. The first model developed to study wind-driven circulation along the northern BC coast was a barotropic, non-tidal, model called the Hecate Strait Model (Hannah et al., 1991; Hannah, 1992). Among the results, it was shown that the wind-driven circulation pattern in HS depends strongly on the earth's rotation. Also, a significant water transport from Hecate Strait was found to be diverted through Chatham Sound on its way into DE. Note that no eddying motions were obtained in Dixon Entrance. Therefore, wind-driven circulation could not explain the presence of eddies in the entrance. Note that a vertically homogeneous water density and velocity were assumed. Such assumption restricts the model application to Hecate Strait only so that model results for other regions are not expected to correspond to observations. In order to deal with the complexity of regional coastline and bathymetry, a model Chapter 2. Review of the Physical Environment and Earlier Models 36 with non-regular grid andfiniteelement techniques was developed (Foreman et al., 1993). It is a barotropic,finiteelement, harmonic tidal model, with realistic bottom topography and coastline, for the entire north coast of British Columbia. It reproduced quite well the elevation amplitudes and phases of the eight largest tidal constituents of the region. However, it was not successful in reproducing the real tidal currents and there was no evidence of the large counterclockwise Rose Spit Eddy. Evidence of a clockwise eddy around Learmonth Bank, with a really complicated flow, appeared. The model has thereafter been improved in order to include a bottom boundary layer (Ballantyne et al., 1996). The tidal elevations and velocities obtained were at least as accurate as those from the previous model. A major difference between the two models is that the Rose Spit Eddy is fully formed in the latter. The period of rotation is however significantly longer than observations; about the same as in the tidal model of Bowman et al. (1992). 2.7.3 Baroclinic numerical models As all the previous models are barotropic, without or with a bottom boundary layer, they were not able to reproduce the significant baroclinic effects that are present in Dixon Entrance. Recently, a baroclinic, diagnostic, finite element model has been developed (Ballantyne et al., 1996). It includes both tidal and baroclinic forcings. In the last few years, another baroclinic model, the Princeton Model (Blumberg and Mellor, 1987) has been applied to the region of interest by P.F. Cummins (Institute of Ocean Sciences, Canada). Thisfinitedifference, time-stepping, model uses a 5x5 km grid which is the same grid used by Hannah et al. (1991). The numerical model used in this thesis is an improved version of the baroclinic, diagnostic, finite element model presented by Ballantyne et al. (1996). A more realistic friction parameterization, a higher grid resolution, the addition of vertical advection and the effect of stratification on friction, as well as the inclusion of wind forcing are among Chapter 2. Review of the Physical Environment and Earlier Models 37 the improvements. The main advantage of this model over the Princeton Model is the easy adaptation of the grid to the complex geography of the region as well as the higher grid resolution. Chapter 3 Description of the Numerical Model and its Physics The governing equations which embody the physics of the problem are the continuity equation and the 3D shallow-water equations of motion for an incompressible and Boussinesq fluid, where turbulent Reynolds stresses are parameterized in terms of an eddy mixing coefficient. These equations, derived in Appendix C, are solved with appropriate boundary conditions in the frequency domain. A wave continuity equation formulation and afiniteelement scheme are used here. The simplest of nonrectangular elements, the linear triangular element, is used. The purpose of this chapter is to describe, in a general way, the model used in the thesis. Note that apart from the parameterization of friction processes, the numerical model used here is the one described in Walters (1992). It is similar to the model used by Lynch et al. (1992), Lynch and Naimie (1993), and Naimie et al. (1996) and employs the numerical methods detailed in Lynch and Werner (1987). Model boundaries, model friction and turbulence, stratification effect, as well as the governing equations with numerical solution procedure are all discussed in this chapter. 3.1 Boundaries of the Model Three physical boundaries need to be considered when studying coastal ocean circulation. Thefirsttwo are the ocean floor boundaries intersect, i.e. the and free land, ocean surface. The curve where these surface is the third boundary. Nevertheless, it is usually de- sired to solve the governing equations for a limited region of interest only. The creation 38 Chapter 3. Description of the Numerical Model and its Physics 39 of an additional boundary, afictitiousone, is therefore necessary. This boundary, commonly referred as open boundary, is the arbitraryfinethat separates the known (model domain) from the unknown (the rest of the ocean). 3.1.1 Ocean floor The bottom boundary layer is generally considered to be made up of three regions; the bed layer (sometimes called viscous sublayer), the constant stress layer (also known as logarithmic layer), and the outer layer. In the ocean, the combination of those three layers are also called the bottom Ekman layer or the planetary boundary layer. Excellent reviews of continental shelf boundary layers are provided by Bowden (1978) and Soulsby (1983). The bed layer, a thin layer adjacent to the sea bed, is typically of a few centimeters thick. The flow in this layer can be classified into smooth, rough, and transitional hydrodynamic roughness regimes. A bed layer with a smooth flow is referred to as a viscous sublayer; a layer in which the bed is sufficiently smooth that the effect of molecular viscosity dominates the dynamics. In general, the total bed shear stress (f ) is the sum of viscous and Reynolds stresses 6 with its x-component, assuming a Boussinesqfluid(Appendix C), given by rl = pJ-^-p (u'w') , t (3.1) t where p» is a reference water mass density, v the kinematic molecular viscosity, u the x-component of water velocity (positive eastward), w the vertical velocity (positive upward), and ' temporalfluctuationsabout a mean (( ) ). In the viscous sublayer, the t viscous stress, i.e. thefirstterm on the right hand side of (3.1), dominates the Reynolds stress, i.e. the second term. The latter can then be neglected and the sublayer velocity profile is obtained by vertically integrating the remaining of equation (3.1). The resulting Chapter 3. Description of the Numerical Model and its Physics 40 velocity is + h) {<)\z v (3.2) where h is the mean water depth and u\ the bottom friction velocity defined as (3.3) Note that the vertical co-ordinate z is defined positive upwards with z = 0 at the mean sea surface. From now on, we will be dealing with mean velocity and the brackets with subscript t will be omitted, unless specified. Typically, the friction velocity is 5% of the velocity value outside of the bottom boundary layer (Soulsby, 1983). If the roughness elements at the bottom boundary exceed about a quarter of the sublayer thickness, which is usually about I2v/u b t (Bowden, 1978), the viscous layer breaks down and the flow is turbulent right down to the boundary; this is the rough regime. In this case, viscosity is unimportant and it is the interaction of the roughness elements with theflowthat determines the velocity and turbulence profiles. The constant stress layer generally extends to a height of a few meters above the bed layer. Neither the details of the bed nor the nature of theflowabove the boundary layer affect its local dynamics. Turbulent mixing dominates and the shear stress is nearly constant and equal to the stress at the boundary. Velocity and turbulence profiles take particularly simple and universal forms. The velocity is constant in direction and increases in magnitude according to a logarithmic profile, that is (z + h) + 6 (3.4) where u is the horizontal water velocity vector, no = 0.4 the von Karman's constant, and £o a length scale associated with the boundary roughness. Finally, there is the outer layer which extends typically to a height of some tens of meters above the constant stress layer. Its velocity and turbulence profiles depend strongly Chapter 3. Description of the Numerical Model and its Physics 41 on the nature of the flow outside the boundary layer, and are therefore not universal. The type of outer layer here considered is one that is produced when a movingfluidflows on a rotating frame, such as on the Earth. The influence of the Earth's rotation causes the current and shear stress to veer progressively clockwise (in the northern hemisphere) with increasing distance from the bed. The velocity vector changes direction as well as magnitude. A bottom boundary layer with this type of outer layer is commonly called a bottom Ekman layer or planetary boundary layer. Generally speaking, the turbulent energy and shear stress decay from a maximum value at or near the bed to zero at the outer edge of the boundary layer. The boundary layer can be defined as the region in which the turbulent energy and shear stress are nonzero. At the outer edge of the boundary layer, the current of the outer layer merges with the geostrophicflow,in which friction is negligible and the pressure gradient is balanced by the Coriolis force. Based on atmospheric boundary layer measurements, the thickness of the bottom boundary layer (A-,) can be expressed as (Davies, 1990) where C\ is a constant and / the Coriolis parameter. C\ of order 0.1 to 0.4 are typical values obtained in atmospheric experiments (Tennekes, 1973). Values of C\ = 0.44 have been found by McPhee and Smith (1976) in the boundary layer beneath drifting ice. More recently, Davies (1990) studied the bottom boundary layer of a coastal region and obtained C\ in the range of 0.3 to 0.6. Note that the values obtained for the ocean are slightly higher than those for the atmosphere. Chapter 3. Description of the Numerical Model and its Physics 3.1.2 42 Free ocean surface Similarly to the existence of a boundary layer at the seabed due to fluid motions, a boundary layer above the water surface exists due to air motions (wind) above the ocean. In response to the wind forcing, a layer, well mixed in mass, momentum, and heat, also develops in the upper ocean. Similarly to (3.5), the thickness of this surface boundary layer (A,) can be expressed as A.=cJ^l, where it* is the friction velocity at the ocean surface. Similarly to (3.6) (3.3), ul is defined as (3.7) with r* being the shear stress at the surface. Using continuity of stress at the air-water interface : f with the wind stress defined as T 3.1.3 W . = f"\ (3.8) Using (3.7) and (3.8), it* can be expressed as Land and open boundaries Along land boundaries, similarly to the formation of the bottom boundary layer, friction between currents and the solid land causes the formation of a lateral boundary layer. However for simplicity, this lateral boundary layer is neglected in the model. In fact as both horizontal viscous and Reynolds stress terms can generally be neglected (see Appendices C and D), this boundary layer can be as well. Chapter 3. Description of the Numerical Model and its Physics 43 Across solid lateral boundaries, i.e. land, depth-averagedflownormal to the boundary is not allowed: (3.10) where rj is the free surface elevation and n a unit vector normal to the boundary. Note that water at a specific depth is however allowed toflowtowards or away from the coast. For example, if surface waterflowstoward (or away) from the coast, the deep water has toflowaway from (or towards) the coast in order to satisfy equation (3.10); this situation corresponds to downwelling (or upwelling) at the coast. The open lateral boundaries are problematic. For now, these boundaries are left aside; they will be treated in the next chapter. 3.2 M o d e l Friction The vertical gradients of Reynolds stresses are important sources or sinks of momentum in the horizontal equations of motion. The Reynolds stresses derive their magnitude from wind stress as well as bottom stress. Frictional processes of a 3D model consequently occur at the free surface, ocean bottom, and in the water column. The model friction at the ocean bottom is computed from the velocity at a given position in the constant stress layer through the use of a quadratic expression that represents the dynamics of the bottom boundary layer. The advantage of this approach is that the dynamics of the bottom boundary layer can be included without having to resolve the entire boundary layer. The friction at the ocean surface is similarly obtained by relating the wind stress to the wind velocity at a specific height above the ocean surface. The internal friction, on another hand, is computed using aflow-dependenteddy viscosity formulation. Chapter 3. Description of the Numerical Model and its Physics 3.2.1 44 Friction at the bottom As a currentflowsover the rough sea bed, frictional stress causes a reduction in velocity, a velocity gradient is thus created, and energy is dissipated through turbulence. Bottom stress is a sink of horizontal momentum, e.g. the bottom boundary layer is where most of the tidal energy is dissipated. By convention, the total bed roughness is partitioned into three parts; the grain roughness, equals to the grain size multiplied by a constant value in the range of 1 to 2.5, the form drag roughness, i.e. ripple or biological roughness, and the roughness caused by sediment motion (Xu and Wright, 1995). The last two are referred to as movable bed roughness. Bottom roughness could then vary in time, but could also vary widely in space from fairly regular features, such as sand ripples, to irregular distributions of objects of sizes ranging from centimeters down to particles of sand and mud. The bottom is a region of no slip. However, bed features can be related to the bed shear stress through a slip condition applied at a distance above the sea bed. One has to be careful because this method does not allow resolution of the physics of the nearbed region, a high-shear area which is known to be important in sediment transport. Nevertheless, as sediment transport study is not the goal of the thesis, a slip condition applied at some distance above the bottom is here considered to adequately represent the physics of interest. The simplest way to relate a slip velocity to the bottom stress is through a quadratic friction law of the form: T = C .(ul+v )*u , B b 2 fP b b (3.11) where Cy is the coefficient of bottom friction, also known as drag or slip coefficient, and u\, the horizontal slip velocity, near the seabed, with components (uf,,Vb). The substitution 45 Chapter 3. Description of the Numerical Model and its Physics of (3.3) into (3.11) gives the relation of the friction velocity and the slip velocity : \<\ = v^J 1*1- ( - ) 3 12 In the constant stress layer, the drag coefficient, Cj, can be viewed as a function of the vertical position; using (3.4), (3.12) becomes ^ ) = { q % j } 2 - < > 313 The value Cj usually used is the Cioo value, where Cioo is the drag coefficient which relates the bed stress to the current 100 cm above the sea bed (this height is in the constant stress layer). From now on, the Cioo value will be used and consequently, the slip currents ui, and Vf, become the currents 1 m (or 100 cm) above the sea bed. The typical magnitude of Cioo is 10 -10 -3 -2 (Csanady, 1976) and in general it is a spatially varying coefficient reflecting bed types and forms. A summary of typical Cioo values for different bottom types can be found in Soulsby (1983). One of the problems in estimating this coefficient is the relative inaccessibility of the boundary layer. The determination of the drag in one location is sometimes considered to be representative of the mean stress over a larger surrounding area of sea bed; the seabed is then assumed to be uniformly rough. As values of Cioo vary through a factor 4 from the smoothest to the roughest substrate (Soulsby, 1983), using a generalized average friction parameterization for an area with different substrates could be inadequate. For example, the residual vorticity generation mechanism by gradients of bottom roughness can be as important as that due to topography (Gross and Werner, 1994). The bottom stress is mainly due to tidal and residual currents. Nevertheless, surface waves can also affect the fluid motion at the seabed. In fact, in water depths less than about half the wavelength of surface waves, surface-wave velocity and pressure fields penetrate down to the seabed. In shallow near-coasted regions, wind waves have then a Chapter 3. Description of the Numerical Model and its Physics 46 significant wave orbital velocity at the seabed producing enhanced levels of turbulence that retards the flow. The presence of surface waves over rough bottoms results in a significant increase in the bottom-friction terms and in additional variability of friction across the shelf. A wave boundary layer, which is much thinner than the current boundary layer seen previously, is formed. The net result is that the current in the region above the wave boundary layer feels a greater resistance than that associated with the physical bottom roughness (Grant and Madsen, 1979). Signell et al. (1990) suggested that wavecurrent interaction could have a significant influence upon theflushingtime of an estuary or semi-enclosed basin. In a numerical model, the appropriate bottom roughness is not just the physical bottom roughness, but is an effective bed roughness that must account for all the energy dissipation processes (e.g. small-scale eddy shedding behind sand waves) that cannot be resolved in the model. The influence of all the unresolved small-scale bed features has to be parameterized. In a 3D model of a near-coastal region, the bottom friction coefficient, and hence the bed stress, has also to be modified to account for the enhanced frictional effects of wind waves at the seabed (Glorioso and Davies, 1995); i.e. wind wave turbulence. These two processes can be parameterized by increasing the drag coefficient; they may be viewed as adding a background drag coefficient. 3.2.2 Friction at the surface At the air-water interface, momentum and energy from the surface wind are transferred to the water. This input generates both low-frequency water motion (currents) and wind waves. For simplicity, only transfer due to the tangential wind stress is considered here; the waves are neglected. Similar to bottom stress, wind stress exerted on the ocean surface is commonly determined by applying the bulk aerodynamic method; the stress magnitude is considered to 47 Chapter 3. Description of the Numerical Model and its Physics be proportional to the square of the wind speed (Geernaert, 1988). The direction of the stress vector is typically assumed to be aligned with the wind. It should however be noted that some observations showed the existence of an angle between wind stress and wind velocity vectors, which depends on height and the atmospheric stratification (Geernaert, 1988; Zemba and Friehe, 1987). Nevertheless, the assumption that wind stress is in the same direction as the wind is made here for simplicity. The bulk parameterization of the wind stress is then defined in a vector form as f = C p (tLl m D + a vl)i* , (3.14) w where Co is the surface drag coefficient, p the air mass density, and u the horizontal a w mean wind velocity at a particular reference height with components (u v ). W) w The drag coefficient is a function of : • sea state (surface roughness), • stability of the air (air stratification over the ocean), • wind speed and fetch, • measurement height of wind speed, and can be written as (Hsu, 1986) C (Z) : D (3.15) ln(^)-i, (Z/L) m where Z is the height above the sea, ZQ the aerodynamic roughness length, L the MoninObukhov stability length, and ip a stratification function: m d(Z/L) 48 Chapter 3. Description of the Numerical Model and its Physics with <f> m defined as Note that in marine environment, the surface waves which affect the roughness length possess a directional distribution. The parameterization of the surface roughness as a scalar is then not exactly right. Observations in the lowest few meters of the atmospheric surface layer showed that the surface wave directional spectrum is important in the directional character of the wind stress (Geernaert et al., 1993). Moreover, Ly (1993) showed that an angle between stress and wind vectors of 10° to 20° would reduce the magnitude of the aerodynamic drag coefficient by 10%-20% for winds less than 30 ms . _1 The effect is even larger for greater values of angle. The ratio {Z/L) is a stability parameter. It is positive in stable and negative in unstable conditions. The stability dependence of Cp can be removed by calculating the coefficient in the equivalent neutral case. Under neutral or adiabatic conditions, L —> oo, Z/L — > 0, <f> (0) m — 1, and thus V'm(O) = 0. The drag coefficient in neutral stratification (CW) is given by Kg CDN{Z) (3.18) and from now on is used. Note that the roughness length, ZQ, is the sum of a roughness height due to shorter surface waves (Z ) and a roughness length for a smooth surface C ( Z ) ; Z — Z + Z, (Smith, 1988). At wind speed > 5 ms , Z is dominant and the _1 B 0 E C drag coefficient increases for increasing winds. At wind speed < 3 ms , Z is dominant _1 A and the drag coefficient increases for decreasing winds. To eliminate its variation with height, the drag coefficient is commonly evaluated at a height of 10 m using the wind speed at 10 m (uio) (Smith, 1988). Using (3.18), it is 49 Chapter 3. Description of the Numerical Model and its Physics Author Smith (1980) Large and Pond (1981) Range of validity (ms *) 6 < u < 22 0.61 + 0.063«io 1.14 4 < u < 10 10 < uio < 26 0.49 + 0.065u Donelan (1982) 0.96 + 0.041u 4 < u < 16 Geernaert et al. (1987) 4 < u < 24 0.577 + 0.085uio Hasselmann et al. (1988) 1.29 0 < uio < 7.5 0.8 + 0.065u 7.5 < u < 50 CW(10 m ) x l 0 3 w 10 lo lo 10 xo lo 10 Table 3.1: Neutral drag coefficients according to various authors (from Blake, 1991). possible to express ZQ as a function of the drag coefficient at 10-m height, such as Z 0 10m = £ X P (\/c (io )) DN Now if CDN (3.19) 1 m is known at 10-m height, the drag coefficient at any height can be obtained using (3.18) and (3.19). The scientific literature suggests a range of neutral drag coefficient as a function of wind speed at 10-m height. Table 3.1 shows various equations which were empirically determined from various oceanic and atmospheric data. The difference between the different values are thought to be due to various sea states under which measurements were taken. Also, observational and model results of several published studies for CDN as a function of Uio are shown infigure1 of Ly (1993). There is relatively little experimental evidence on the behaviour of the drag coefficient at wind speeds below 3 ms . Further-1 more, estimation of wind stress at wind speeds above about 25 ms remains speculative, -1 since few if any direct measurements are available. Following the suggestion of Large and Pond (1981), the neutral drag coefficient is taken here as C (10 DN m) x 10 3 = 1.2 = 0.49 + 0.065uio (Ams' <u < 1 10 (11 ms -1 <u Urns' ), 1 < 25 ma" ). (3.20) 1 10 Chapter 3. Description of the Numerical Model and its Physics 50 However, since most of the wind speed values observed are less than 11 m s , a simplified - 1 version of (3.20) is actually used in the thesis. CDN{10 rn) is considered as a constant and set to 1.2 x 10~ for all wind speed values. 3 3.2.3 Friction in the water column Across a horizontal element of area perpendicular to the z-direction, the vertical transport of horizontal momentum downwards is equivalent to the total horizontal shear stress by which the water above the area acts on the water below; it is given by equation (3.1). The bottom friction is parameterized here at a height of 1 meter above the seabed, i.e. in the constant-stress layer. Above this height, the contributions to the vertical transport of momentum by molecular processes is negligibly small compared with the contributions of turbulent eddies. The shear stress can therefore be related to the gradient of mean velocity by an eddy viscosity assumption of the form, for x-component, r\ = -p.{u'w') = P.A ^, t (3.21) V where A is called the vertical eddy viscosity coefficient. Unlike molecular viscosity, A v v varies with position and with the mean flow or the level of turbulence in a way which has to be assumed. Omitting the brackets with subscript t, (3.21) becomes rl = P.A —, (3.22) V where u is the mean velocity. The simplest assumption is to take A v as constant in space and time. However, constant eddy viscosities over depth are incapable of representing the bottom boundary layer in a physically realistic manner. Also, its use in a region of significant water depth variation, and changes in current, is not physically realistic (Davies, 1986). We are therefore looking for a better approximation than constant eddy viscosity. A is actually v Chapter 3. Description of the Numerical Model and its Physics 51 the momentum exchange coefficient in the vertical and can be defined by an eddy length scale times a RMS eddy velocity. Near the sea bed, this eddy length scale is proportional to K (Z + h) and the eddy velocity scale to the friction velocity. Following the procedure 0 by Grenier et al. (1995), Davies (1993), and Aldridge and Davies (1993), A is, in the v bottom boundary layer (A ) , defined as b A = K l {z)\ut\, (3.23) b v 0 b where h{z) is a mixing length still to be determined. One advantage of (3.23) is that •u^ is derived directly from the computed bottom stress, rather than from the velocity profile. Using (3.12), equation (3.23) can be re-written as A = K (C )h (z)\u \. (3.24) b v 0 100 b b An approximate model for steady flow in a bottom boundary layer can be constructed using an eddy viscosity that increases linearly near the sea bed and is constant throughout the rest of the boundary layer (Smith-Long, 1976). In Grenier et al. (1995), h(z) is taken as increasing linearly with height above the sea bed over the lower 20% and remaining constant over the upper 80% of the bottom boundary layer. This assumption was based on measurements made in tidal currents by Bowden et al. (1959). In the model used here, the mixing length is defined in the same way as in Grenier et al. (1995), except that h(z) is taken as decreasing linearly with height over the top 20% of the boundary layer. In this way, l (z) reaches zero at the boundary limit. The mixing length in the bottom b boundary layer is therefore expressed as l b = z+ h - h < z < -h + 0.2A , 6 h = 0.2A - h + 0.2A < z < -h + 0.8A , l = A - (z + h) -h + 0.8A < z < -h + A , l = 0 - h + A < z. 6 b b 6 6 6 b 6 6 (3.25) Chapter 3. Description of the Numerical Model and its Physics 52 The bottom boundary layer thickness (A-,) can be re-expressed, using (3.5) and (3.12), as A & = C i ( ^ o o ) i . | | ( 3 2 6 ) Afc then evolves with theflowfieldand its magnitude changes with horizontal position. In coastal studies, it is often the case that the water depth is less than the thickness that the boundary layer would otherwise attain. The growth of the boundary layer is thus restricted and the whole depth of water is turbulent. Therefore, if the computed bottom boundary layer exceeds the total water depth (H = h + n), A-, is depth-limited and set to H. Similarly to (3.23) for the bottom boundary layer, an eddy viscosity formulation may be obtained in the surface boundary layer: A' V = (3.27) K 1,{Z)\U:\. 0 Using (3.9) and (3.14), (3.27) now becomes A' V = K 0 (JCDN^J 2 l.{z)\u \, w (3.28) with l (z), similarly to (3.25), given by a l a = 0 z < - A , + rj, I. = A + (z - n) - A + n < z < -0.8A, + n, l = 0.2A, - 0.8A, + n < z < -0.2A, + rj, l = - z + 77 - 0.2A, + 7/ < z < n. a a a a (3.29) Note here again that if the total water depth is shallower than the surface boundary layer, A , is depth-limited and set to H. For typical nonstorm conditions in mid- and outer-shelf regions, the surface and bottom boundary layers are well separated by a non-frictional core region. The vertical Chapter 3. Description of the Numerical Model and its Physics 53 eddy viscosity coefficient is then equal to (3.24) in the bottom and to (3.28) in the surface Ekman layers. However during storms in shallow waters, or in waters with high tidal currents, the two boundary layers can merge. By summation of the turbulent kinetic energy due to the two processes, the total eddy viscosity coefficient (A°) might be expressed by (3.30) Nevertheless, for model simplification, A° is taken here simply as -40 _ Ab i AS (3.31) Note that tidal turbulence at depth is important in diffusing the wind's energy out of the surface boundary layer. In regions of strong tidal currents, tidally induced turbulence at depth is large and the wind's momentum diffuses very rapidly from the surface layer. Actually, the wind-driven currents are found to be smaller in regions of high tidal current than in regions of low tidal current. Davies (1985), and Davies and Lawrence (1994) showed that in shallow water where tidal turbulence could reach the surface, the winddriven currents were weaker than in deep water. The difficulty of choosing an eddy viscosity profile and formulation can to a certain extent be removed by using a higher-order turbulence closure, that is a turbulence energy closure scheme. However, the advantage of higher-order closure models, based on prognostic equations for turbulent energy and its dissipation, at present remains unestablished in shallow-sea modelling (Aldridge and Davies, 1993). Davies (1991) found that the uncertainty in some of the parameters associated with the mixing-length formulation were comparable with the uncertainty in the knowledge of a simpleflow-dependenteddy viscosity model. Also, Xing and Davies (1995) found that simpleflow-dependenteddy viscosity models do as well as turbulence kinetic energy models. Note that three friction parameters need to be introduced here into the model; these are Ciooj CDN, and C \ . Chapter 3. Description of the Numerical Model and its Physics 3.3 54 Effect of Stratification The water column on the continental shelf generally exhibits vertical gradients in density associated with gradients in temperature, salinity, and (in special cases) suspended sediments. When afluidis stably stratified, shear production of turbulence is partially balanced by work done against buoyancy forces: the level of turbulence is reduced relative to the unstratified level. The boundary layer is reduced in thickness and the vertical extent of the turbulence is constrained to scales smaller than the neutral Ekman layer thickness. Thereby vertical mixing is reduced and so is the vertical eddy viscosity. The net result is that the transport of momentum across the sharp density interface is inhibited to some extent and so is the diffusion of properties such as temperature. When stratification occurs, an internal adjustment will also occur in the pressurefield;baroclinic motions may be generated: internal tides, density currents, etc. Such motions can also affect the vertical eddy viscosity because they change the vertical current structure and thus the vertical mixing in the water column. The stratification of surface waters depends not only on the conditions at the surface but also on the energy generated near the seabed. In shallower water on the continental shelf, turbulence generated in the bottom boundary layer by tidal currents is an important factor in determining whether a seasonal thermocline develops. If sufficient turbulent energy is transported upwards from the bottom layer it may prevent the formation of a thermocline and maintain vertical mixing through the water column (Simpson and Hunter, 1974). In order for mixing to take place throughout the water column it is necessary not only that there be enough turbulent energy available to overcome the potential energy associated with the stratification, but also that the energy, which is generated near the seabed, be capable of reaching the water surface. A number of stability parameters exist. There is the Brunt- Vdisdld frequency (N), 55 Chapter 3. Description of the Numerical Model and its Physics also known as stability or buoyancy frequency, which measures the static stability. For an incompressible Boussinesq fluid (Appendix C), it is defined as N 2 = _g_0fo p* oz ( 3 3 2 ) with z positive upwards, g the gravitational acceleration (g =9.8 ms ), and po(z) a -2 mean water density. If N > 0, the water is stable and a parcel displaced a short 2 distance vertically will tend to return to its original position. Because it has inertia, it will tend to overshoot its original position and then to oscillate about it, hence the stability of the water is related to the properties of internal waves. If N = 0, the water 2 is neutrally stable and a displaced parcel will tend to remain in its displaced position. If N < 0, the water will be unstable and a parcel which is displaced will tend to continue 2 its displacement, i.e. overturn of the water should occur. Another stability parameter is the dimensionless Richardson Number (Ri), sometimes called the gradient Richardson Number, which measures the relative importance of mechanical mixing and density stratification. It is given by (3.33) If Ri < 0, density variations enhance turbulence; if Ri > 0 they tend to reduce it. The stability is commonly defined as near-neutral if \Ri\ < 0.03. In order to calculate Ri, simultaneous velocity and density profile measurements are required. A third parameter is the dimensionless flux Richardson Number (Rf) which expresses the ratio of energy consumed by buoyancy forces to stress production of turbulent kinetic energy. This number involvesfluctuatingquantities that are difficult to measure in practice, so the gradient Richardson Number is generally used in applications. Simple formulas relating the eddy viscosity coefficient to the Richardson number are Chapter 3. Description of the Numerical Model and its Physics 56 found in the literature. Munk and Anderson (1948) define A = F(Ri)A° , v (3.34) v where A„ is the vertical eddy viscosity for neutral stability defined in (3.31). F(Ri) is a reduction factor crudely representing the influence of stratification in terms of a time-independent overall Richardson number. The general form for F(Ri) is where X\ and X2 are adjustable positive parameters. Ruddick et al. (1995) found that this type of simple Richardson number dependence for diffusion coefficients can give results which aire qualitatively similar to those of more complicated models. However since less of the physics of turbulence is included in the simpler models, more calibration is required. 3.4 Governing Equations of the Model In this thesis, the phenomenon of interest is the mean circulation. The solution of interest is then expected to change at a slow enough rate for the acceleration terms to be negligible; that is a steady-state solution. Indeed, in the frequency domain, a solution for the null frequency (a; = 0) is sought; a diagnostic model is considered. The set of governing equations for this steady-state study, considering weak nonlinearity and neglecting for now the stratification effect on turbulence, is (see Appendix D for details on their derivation) : V •u, +^ h V h fx Ures + g^res re • (H ^ ) m res oz = V h = 0, (3.36) •Q , (3.37) res " ^ f ^ ^ f ) = ^s + fres + , (3.38) Chapter 3. Description of the Numerical Model and its Physics / x u , + gV <; + re h = *™ + Rre. + T rea rea 57 - (3.39) Um rim with boundary conditions : + f„. = H % m ea {z = 17), (3.40) ~h ), (3.41) Oz &1L -* A l ^ - -*• + Tres = C 1 0 0 l U r e T + s (Z = r e i m OZ 1 N Wre, = {U . • Vfcfc. .) + - £ re 4 Wre. = - Ures • ft = («P ' ^ - p ) («re« (* = »?)> V e ( 3 4 2 ) p=-JV * Vfc/l ) (z = - / l ) , m m ((z,y) 0 = 'and b o u n d a r y ) , (3.43) (3.44) first-order vertical viscosity : Al = K (C7 O)"/6(Z)7 + ° K 0 10 { H ^ r e . \ f /.(*), (3.45) nonlinear forcing terms : Qre. = " J £ %3-p, (3-46) P^O T J- - - re* — 11(^(^1 + ^ ) , (3.47) p=-N p/o and high-order friction terms P^O f ^ ^ K ' O ^ W For any variable X, X Tes (3-49) p=-JV P^O is the residual component of X while X is the tidal component p with frequency w of X and X-p the complex conjugation of X . Note that N defines the p p Chapter 3. Description of the Numerical Model and its Physics 58 number of tidal constituents considered in the study and X indicates the vertical average of X. I should point out here that the bottom of the model is set to 100 cm above the sea floor. At this defined bottom, thefluidis constrained to follow the topography and the bottom stress is a sink of horizontal momentum, iib, named bottom velocity, is the velocity at the model bottom, c the amplitude of the free surface elevation, h the m mean water depth of the model, defined as h = h— 1 m, and H the total water depth m + c ). / is the Coriolis parameter, H ^ , (H = /i Rret, the baroclinic pressure gradient vector (see m m re5 of the bottom speed (see 3.5 m m re the atmospheric forcing (see D.13), and 7 D.10), the time-independent part D.21). Numerical Solution Procedure The model requires known forcings which are provided here by tides, wind, baroclinic pressure gradients, and river runoff. The model solution is actually a linear combination of fluid velocities due to each of those driving mechanisms under specific open boundary conditions. Note that atmospheric pressure gradients are assumed negligible over the scale of the region of interest. 3.5.1 Solving the governing equations The computational task is to solve the continuity and momentum equations subject to external forcings as well as forcing from the nonlinear tidal terms. Tidal residuals are exclusively generated by nonlinear processes and thus inclusion of advection effects is important. The interaction term of the residual flow with itself is however assumed negligible and neglected here (weak nonlinearity). Prior to being discretized, the continuity and momentum equations are combined into a Generalized Wave Continuity Equation (GWCE) which has been shown to have superior 59 Chapter 3. Description of the Numerical Model and its Physics numerical properties to a primitive continuity equation when afiniteelement method is used in space (Lynch and Gray, 1979). The GWCE is solved together with the primitive momentum equations using a Galerkinfiniteelement method with nodal quadrature for the elevation of inner products on linear triangles in space. Extensive details of the numerical solution strategy can be found in Appendix D. Here is a summary of the main sequential steps: . —* 1. Surrogate variables, E = (E ± res —¥ • x ± iE Tea • y)/2, are defined in order to uncouple the momentum equations (3.38) and (3.39) (see Section D.6). 2. Equation (3.38), a one-dimensional (ID) vertical diffusion equation, is solved with boundary conditions (3.40) and (3.41) for u (z) in terms of two partial sores lutions and the unknown c ,. Pe 3. The bottom velocity (uj,.,), obtained during the preceding step, is substituted into equation (3.39) which is then substituted into the continuity equation (3.37). This creates the GWCE (equation (D.77)), a Helmholtz equation for c . Pe4 4. The resulting Helmholtz equation is solved with boundary conditions (3.44) for c . This solution provides the barotropic pressure response which accompaPe4 nies the imposed density field, tides, wind, and open-water barotropic boundary conditions (including river runoff). 5. The horizontal velocities, u ea, are reconstructed using the two partial solutions T obtained during step 2 and the known c ,. Pe 6. Finally, the 3D continuity equation (3.36) is solved with boundary conditions (3.42) and (3.43) for the vertical velocity w res using the known u . As the rea continuity equation is an ordinaryfirstorder differential equation, it can satisfy Chapter 3. Description of the Numerical Model and its Physics 60 only one boundary condition. The problem is that we have two boundary conditions available here. In order to solve this overdetermined system of equations, the least squares method presented in Muccino et al. (1994), and also in Muccino et al. (1996), is used. In brief, the solution of the governing 3D equations simplifies to a ID vertical diffusion equation for the horizontal velocities (horizontal momentum equation), followed by a twodimensional (2D) horizontal Helmholtz equation for the elevation (GWCE), and finally a 3D solution of the continuity equation for the vertical velocity. These three equations are solved using thefiniteelement method detailed in Appendix E. 3.5.2 Relaxation and convergence test Since bottom friction and turbulence in water column (eddy viscosity profile) depend on the bottom velocity, which is not known a priori, an iterative scheme is required. To begin the iteration, an initial guess of the bottom velocity is made based on tidal currents. This estimate is used to determine the bottom stress and eddy viscosity profile used in the calculation of the new bottom velocity, and then in the calculation of surface elevation and horizontal velocities. Relaxation is used to counteract any instabilities which may occur from the nonlinear terms. The amount of relaxation is decreased as the friction parameters and steady velocities converge. After each iteration of the solver, the solution is updated using : $ = $,+ - r r ( $ r - $ + ) , i where = updated solution at node i , $7= solution at node i prior to the latest iteration , = solution at node i computed during the latest iteration , r = relaxation parameter . (3.50) 61 Chapter 3. Description of the Numerical Model and its Physics $ i represents either the bottom or depth-averaged velocity. The relaxation parameter is set to 0.2. After each run, the surface elevation and depth-averaged velocities are compared to those just prior to the run using the following : NN El*r-*?l A* = '— v e NN *L. , (3.51) = max(\*7-*t\), (- ) 3 52 where NN is the total number of nodes in the horizontal. In this thesis, convergence is met when : A * < 0.02 cm s e A* ve < 0.02 cm A^ < 0.1 cm s a!B A* .< 0.1 cm aa _1 for $ = depth-averaged velocity , for $ = surface elevation , _1 for $ = depth-averaged velocity , for $ = surface elevation are all satisfied. This usually takes between 10 and 30 iterations. Chapter 4 Modelling the Dixon Entrance System Once the model's ability to solve the governing equations has been thoroughly examined, one may proceed with confidence to the substantive issue of modelling a real system. This chapter describes the spatial discretization of the domain, the choice of open boundary conditions and model parameters, as well as the determination of model forcings. Only forcings from tides, baroclinic pressure gradients, wind, and river runoff are part of this thesis. Model results as well as evaluation of model sensitivity to variations of grid resolution, open boundaries, and friction parameters are left aside for now; they will be the subject of the following two chapters. 4.1 Spatial Discretization A physical coordinate space 3D Finite Element Mesh is used. It is created by projecting a 2D (x,y) mesh, composed of linear triangular elements, vertically from the sea surface to the sea floor. The vertical mesh, consisting of linear elements, is constant under a specific horizontal node but varies with horizontal position. The number of vertical nodes is constant through the domain. This 3D mesh is next transformed from the physical (x,y,z) to a mapped (x,y,(r) coordinate space (see Fig. 4.1). Wetting and drying elements are not allowed here. In order to prevent any element to dry, a minimum bathymetric depth of 5 m is imposed. The bathymetry data used come from different sources. First, a topographic database (~lxl km) available from Energy, Mines and Resources Canada (EMR) was used. 62 Chapter 4. Modelling the Dixon Entrance (a) 63 System (b) Figure 4.1: 3D mesh in (a) physical (x,y, z) and (b) mapped (x,y,cr) coordinate space. The data were smoothed and interpolated on the corresponding triangular grid by L . E . Cuypers (Victoria, B C ) . As only the Canadian waters are included into this database, the bathymetry for American waters was obtained from National Oceanic Atmospheric Administration ( N O A A ) bathymetric charts. In Canadian waters, gaps were also present. In that case, Canadian Hydrographic Service bathymetric charts were used. I also compared the bathymetry obtained from the database with charts and corrected it when inaccurate. The database from which the coastline data were extracted comes from the National Geophysical Data Center (NGDC), a division of N O A A . 4.1.1 Horizontal discretization The horizontal domain is discretized using a finite element mesh of linear triangles generated with the T R I G R I D software described by Henry and Walters (1993). This grid permits a good fit of the irregular coastline and bathymetry characterizing the study Chapter 1 3 4 ' 4. the Dixon 133 W w ~ 134 W Modelling r 133 w Entrance 64 System 132 W , 132 W 131 W , 131 W 130 W — 1 130 W Figure 4.2: Coarse G r i d - Horizontal triangular mesh composed of 1742 elements and 1038 nodes. Numbers represent the different open boundaries. area. Elements of an arbitrary size, shape, and orientation are used. Smaller elements are placed in regions where either the bathymetry, surface elevation, or velocity field is changing rapidly and higher accuracy is desired. For this research work, four different grids were used. The first one, called Coarse Grid, is characterized by a resolution ranging from 16 km to 1.2 km with 76% of the grid lying in the 4 km and higher range (Fig. 4.2). The second grid, named Fine Grid, is a refined version of the Coarse Grid. It has a resolution extending from 8 km to 0.8 km with 86% of the grid situated in the 4 km and lower range (Fig. 4.3). This finer grid is used in order to ascertain that grid resolution is fine enough so that the introduction of any Chapter 4. Modelling the Dixon Entrance System 65 Figure 4.3: F i n e Grid-Horizontal triangular mesh composed of 5432 elements and 2982 nodes. Numbers represent the different open boundaries. Chapter 4. Modelling the Dixon Entrance System 66 major numerical errors in the solution is avoided. Hannah and Wright (1995) effectively found errors occurring at the coast and in regions with abrupt changes in bottom slope when using a similar numerical model. Thefirstmodel weakness was related to the fact that the estimate of pressure gradient at the coast is essentially one-sided and thus only first order accurate. The second weakness occurred in the pressure gradient estimate due to the interpolation of poorly resolved length scales in the pressure field. It was suggested that these two errors might be corrected by using either a higher order scheme or by adding extra horizontal resolution. In order to evaluate the impact of the open boundaries on the numerical solution, an extension of the previous grid was created (Fig. 4.4). It is called Extended Grid and has the same spatial resolution as the Fine Grid. The only differences between these two grids are open boundaries 4, 5, and 6 which have been pushed away in the Extended Grid. All the previous grids are in reality subdomains of larger grids which are used to model tides. Note however that tidal modelling is not part of this thesis. The Coarse Grid is a subdomain of the grid used in Foreman et al. (1993). The Fine and Extended grids are both subdomains of the Large Grid shown in Fig. 4.5. The latter is also used in the modelling process of tidal residuals and wind-driven currents. 4.1.2 Vertical discretization Using the ordinary z-level coordinate system (Fig. 4.6) has certain disadvantages in the vicinity of large bathymetric irregularities. For example, the deep region in Fig. 4.6 is discretized with as many as seven elements while the shallow region is discretized with only one. Note that the bottom of the modelled domain is displaced to a small distance above the sea floor, taken as 1 m. It is desirable to introduce a terrain-following acoordinate system in which the water column is divided into the same number of grid Chapter 4. Modelling 134 W 134W the Dixon 133 W 133W Entrance 67 System 132 W 132W 131 W 131 W 130 W 130W Figure 4.4: E x t e n d e d Grid-Horizontal triangular mesh composed of 6128 elements and 3350 nodes. Numbers represent the different open boundaries. Chapter 4. ModeUing the Dixon Entrance System 68 Figure 4.5: Large Grid-Horizontal triangular mesh composed of 13,133 elements and 7575 nodes. 69 Chapter 4. Modelling the Dixon Entrance System Figure 4.6: Vertical grid with a z-level coordinate system in physical (x,y,z) coordinate space. cells independently of horizontal position (Phillips, 1957). This system can thereafter be transformed in a map (x,y,a) coordinate space where the surface and bottom of the model are set to a = 0 and a = — 1 respectively. The interval rj > z > — h , which varies m with horizontal position, can then be rewritten as 0 > z* > — H m = —(h + 77), with m z* = z — 77, and thus converted into the constant interval 0 > cr > —1. Two vertical grids, characterized by 10 equally-spaced linear elements (11 nodes) and 20 equally-spaced linear elements (21 nodes), were used (Fig. 4.7). It is matter of interest to note that grid resolution (Az) varies with water depth; high vertical resolution characterizing shallow regions and low resolution deep areas. The hydrodynamic equations (3.36-3.49), written in Cartesian coordinates, should also be converted to sigma coordinates. Define the Cartesian and sigma coordinate systems as (x,y, z) and (a;*,y*, o~) respectively. The transformation from one to the other Chapter 4. ModeUing the Dixon Entrance System 70 Figure 4.7: Vertical linear meshes, in physical (x,y,z) coordinate space, with (a) 10 elements and 11 nodes and (b) 20 elements and 21 nodes. system is given by: z-y{x,y) x =x, y =y, o- = -rr-,—r-- (4.1) Hm{x,y) The vertical derivative in sigma coordinates is then given by 1_ _ 1 1 (4.2) d~z ~ H~Zd~a' and the horizontal derivatives by 1 1 dx dx* 1 1 dy dy* 1 dn H m dx 1 drj Hm dy + 0- dH„ dx J_ dH m d^' dy (4.3) can be approximated by 1 1 a dh 1 m rsj dx ~ : dx* h dx da- 1 m 6V (4.3) Chapter 4. Modelling the Dixon Entrance System 1 1 dy dy* — ~ o- h m dh 71 m 1 dy do- . (4.4) Note that the first term on the right hand side of (4.4) involves the horizontal gradient along a constant <r-surface while the second term involves the gradient of the bottom topography. When calculating the pressure gradient in sigma coordinates near steep topography, these two terms become large, comparable in magnitude, and opposite in sign. In such a case, a small error in computing either term near steep topography can result in a large error in the total pressure gradient force (Haney, 1991). This is however not a problem here since the baroclinic pressure gradients are calculated directly from observations and not through model calculations. Note that the model calculates the barotropic pressure gradient (V^r?) but in this case the second term on the right hand side of (4.4) is null (g = 0). 4.2 Open Boundaries Open boundaries are fictitious boundaries created by the fact that only the water circulation of the Dixon Entrance System is modelled. The only way to eliminate these boundaries would be to consider the entire earth as the model domain. This process is however too challenging for this thesis. The role of open boundaries is to link the physics of the study area to the physics of the rest of the world. If no open boundaries were set up around the model domain, then no water would be allowed toflowin and out of the study region. In that case, only the currents generated inside the model domain would show up. In addition, unphysical currents might appear in regions, close to boundaries, where water is expected toflowin or out of the domain. The Dixon Entrance System is characterized by numerous open boundaries (see Figs 4.2 to 4.4). Any water transport through these boundaries needs to be specified. The maun problem is to choose the right boundary conditions. In the model used here, either Chapter 4. Modelling the Dixon Entrance System 72 depth-averaged currents (direction and/or intensity) or surface elevation can be specified at boundaries. The values can be obtained either from observations or from modelling a larger domain that includes the study area. Based on observations (see Chapter 2), the direction of depth-averaged current is set to be normal to all open boundaries except for boundary 5 where it is set to be tangential. No specific value for the depth-averaged water speed is imposed on open boundaries. The only exception is at river mouths (boundaries 11 and 13) when the river runoff forcing is present. In this case, the speed is directly obtained from the riverflowratedischarges (see Section 4.6). As at least one specified elevation is required to permit a solution to the GWCE (wave equation with Neumann boundary conditions), the elevation of the node located at the southwest corner of the grid is set to zero. Note that only the relative elevations are matter of interest in the model calculations and therefore the particular node and its specified value are both arbitrary. When running the model on the Large Grid, we can model the currents of the DE System without having to bother about open boundary conditions. On this grid, the open boundaries are pushed so far away from the region of interest that the influence of the boundaries on the solution is negligible. In this case, the surface elevation is then simply set to zero at all open boundary nodes. It should however be mentioned that it is not possible, due to lack of observations, to model the density-driven currents on the Large grid. 4.3 4.3.1 Model Parameters Friction parameters Energy enters the area through different processes: wind forcing at the sea surface, tidal forcing at ocean boundaries, and density-driven forcing at river mouths. There 73 Chapter 4. Modelling the Dixon Entrance System axe two dissipation mechanisms for removing this energy: bottom friction, the primary mechanism, and frictional processes within the water column. The representation and parameterization of dissipation mechanisms remain a delicate aspect of modelling. In this thesis, frictional processes are computed using a wind stress-velocity relationship above the sea surface (3.14), a slip-condition at 1 m above the seabed (3.11), and in the interior, a combination (3.31) of a flow-dependent eddy viscosity for the bottom boundary layer (3.24) and a wind-dependent eddy viscosity for the upper boundary layer (3.28). The primary means of calibration for most hydrodynamic models is the adjustment of friction parameters, taken here as CDN, CWO, and C\. The former, CDN, has already been discussed in Section 3.2; its value is fixed to 1.2 x 10~ at 10-m height. The other 3 two parameters are still to be determined. Little is known about the nature of sea floor sediments (bed types and forms) and of flow turbulence in the region under investigation (see Section 2.2). The choice of bottom and interior friction parameters can then hardly be based on observations. As in Ballantyne et al. (1996), the bottom friction coefficient (Cj = Cioo) is set to 0.01. This value was used to model tides in a region that included the Dixon Entrance System. I should however point out that the representation of the frictional processes within the water column used by Ballantyne et al. (1996) is different from that detailed in Chapter 3. The difference lies in the definition of the bottom boundary layer thickness (A&) and the mixing length in this layer (lf,(z)). It should also be kept in mind that using a spatially constant C o value over the entire model domain does not reflect the bottom variations 10 of the region of interest. For now, it is considered only as a first-order approximation. The third friction parameter to fix is C\. Note that the thickness of boundary layers ((3.5) and (3.6)), and thus the mixing length in these layers ((3.25) and (3.29)), depend on this parameter. C\ is set here to 0.3 which is a value falling into the range of values obtained from different studies (see Section 3.1). 74 Chapter 4. Modelling the Dixon Entrance System As in Naimie et al. (1994), a background vertical eddy viscosity, A b a = 0.002 m s , is added. The vertical eddy viscosity under neutral stability (3.31) is now 2 - 1 expressed as A° = A b v + A> -r A J> (4.5) b v v and (3.45) replaced by Al = Ko(C 1 0 o ) ^ ( z ) 7 + "0 ( # n . | * r e - . | ) ' + A*. (4.6) The addition of this constant term can be viewed as simulating the contribution of unmodeled currents. 4.3.2 Other parameters Among the other parameters that one has to provide the model with is the Coriolis parameter / . It is defined as / = 2VLsin<f> with = 7.292 x 1 0 - 5 s - 1 and <p the latitude (see Appendix C ) . The latitude chosen for this study is 54.3° N . It approximatively corresponds to the value at the center of the modelled domain: the Dixon Entrance System. In order to calculate the baroclinic pressure gradient (R) and the wind stress (H if>), a reference density for sea water (p„) and for air (p ) are needed (see (D.13) m a and (D.10)). We are using the values suggested by Gill (1982) : p, — 1025 kg m - 3 and = 1.2 kg m - . 3 P a A summary of all model parameters used is presented in Table 4.1. Note that the stratification parameters, Xi and X , are missing. 2 parameters in the next chapter. We will come back to these two Remember that for now the stratification effect on friction is neglected and therefore not incorporated into the governing equations. 75 Chapter 4. Modelling the Dixon Entrance System Value Parameter Von Karman's constant Gravitational constant Angular velocity of Earth Latitude Coriolis parameter Reference density of sea water Density of air Neutral surface friction coefficient Model bottom friction coefficient Coefficient for bottom and surface boundary layers Background eddy viscosity Ko 9 n f P* Pa CDN ClOO Ci 0.4 9.8 m s 7.292 x I O s54.3° N 1.184 x 10" s" 1025 kg m " 1.2 kg m " 1.2 x 10" 0.01 0.3 0.002 m s " - 2 - 5 1 4 1 3 3 3 2 1 Table 4.1: Model parameters. 4.4 Tides The inclusion of tidal constituents in a steady-state model is important. First of all, the neglect of tidal constituents would mean that both the eddy viscosity and bed stress computed with the model would be smaller than those found in nature, since both terms depend upon the total current and not just the steady-state flow. Corrections could simply be made by increasing the viscosity and bottom friction but the spatial variation of these terms would be hard to determine. The second reason is that important tidal residual currents are generated by nonlinear processes. If the tidal advection effects were neglected, these processes would not be modelled. In this perspective, tidal constituents are permanently included in the friction terms of the model but not in the forcing terms. The nonlinear terms need to be present in order to have tidal forcing. Chapter 4. ModeUing the Dixon Entrance System Tidal Constituent M 2 s o N 2 1 2 Pi K 2 Qi Frequency (cph) 0.0805 0.0418 0.0833 0.0387 0.0790 0.0416 0.0836 0.0372 76 Period (hours) 12.42 23.92 12.00 25.84 12.66 24.04 11.96 26.88 Table 4.2: Tidal constituents, in order of importance, included in the model. 4.4.1 Tidal constituents In coastal waters, tides are generally classified as astronomical tides, compound tides, and overtides. Astronomical tides result from the gravitational forces exerted by the sun and the moon. Compound tides and overtides, sometimes called shallow-water tides, arise from the nonlinear interactions between astronomical constituents. More specifically, compound tides result from interactions between two or more constituents of different frequencies, whereas overtides result from interactions of a single constituent with itself or its overtides. For simplicity, only the first eight most important tidal components are considered in the model; they all are astronomical tides (see Table 4.2). These eight constituents account for 89% of the tidal height range in the daily and half-daily frequency bands at Prince Rupert (Foreman et al., 1993). The M constituent is produced by the 2 moon orbiting around the earth and S by the earth orbiting around the sun. It takes 2 24 lunar or solar hours for the earth to expose the same point to the moon or the sun and the bulge of tidal forces pass twice through this point during the time period. The inclination of the earth's axis and the moon's orbit introduces a complication from which the Ki, Oi, Pi, and K constituents arise. A further complication is introduced by the 2 Chapter 4. Modelling the Dixon Entrance System 77 variable distance of the moon and the sun from the earth. This produces the N and Qi 2 constituents. Note that no compound tides or overtides are considered in this study. However, it should be pointed out that a nonlinear interaction between K\ and 0\ contributes to the M tide; similarly, K\ and P contribute to S , and K\ and Qi to N . These nonlinear 2 x 2 2 interactions have been included in the modelling of astronomical tides. 4.4.2 Calculation of tidal constituents The modelling of tidal constituents is not part of this thesis. Tidal elevations and horizontal velocities were provided by M.G.G. Foreman (Institute of Ocean Sciences, Canada). Refer to Foreman et al. (1993) and Ballantyne et al. (1996) for details on the tidal model. Note that the elevations and horizontal velocities were calculated here using a modified version of the tidal model presented in Ballantyne et al. (1996). For consistency, the parameterization of frictional processes within the water column was replaced by that presented in Chapter 3. In this way, both the tidal and diagnostic models are using the same physics. The vertical tidal velocities were calculated in a similar way to the calculation of the vertical residual velocity presented in Chapter 3. For a specific tidal constituent, the 3D continuity equation (D.22) and boundary conditions (D.33) and (D.34) are respectively given by dw = 0 dz (4.7) p W p = -IWpSp -Up • V h, 'm (4.8) {z = rj) and w = p h -h ). m (4.9) Chapter 4. Modelling the Dixon Entrance System 134W 133W 134W 133W 132W 132W 78 131 W 130W 131 W 130W Figure 4.8: Depth-averaged nonlinear tidal forcing term (T ,(x,y,z)) equation. Units are cm s and values were calculated on the Large re - 2 in the momentum Grid. The vertical tidal velocity w is then obtained by solving (4.7) with (4.8) and (4.9) using p the known u and q . Note that for simplicity the non-linear term is neglected here. p 4.4.3 p Tidal forcing Tidal forcing enters in the model via two nonlinear terms. The first one, given by (3.47), constitutes the advection of tidal momentum by the tide. Its depth-averaged value is shown in Fig. 4.8. The second nonlinear term, given by (3.46), corresponds to 79 Chapter 4. Modelling the Dixon Entrance System 134W 134W 133W 133W 132W 132W Figure 4.9: Nonlinear tidal forcing term (Q e») m s and values were calculated on the Large r 2 _1 131 W 131 W 130W 130W in the continuity equation. Units are Grid. the wave transport (see equation (D.77)) (Walters, 1992). This is equivalent, in the onedimensional case, to the Stokes drift. This tidal forcing term is shown in Fig. 4.9. Notice that in interest of clarity, the computed forcing terms were interpolated onto a 4x4 km grid and then plotted. This interpolation was applied to all forcing terms presented in this chapter. Chapter 4. Modelling the Dixon Entrance System 4.5 80 Baroclinic Forcing The different baroclinic pressure gradients, used to drive the numerical model, were calculated from density fields characterizing the study area at different time periods. These density fields were obtained, with the help of statistical interpolation, from CTD observations taken during different cruises in the region of interest (see Appendix B for a list of cruises). Note that density data have not been directly measured but rather calculated, with a well establish formula (Gill, 1982), from temperature, salinity, and pressure measurements. Before using them, the density data were checked in order to eliminate any clearly erroneous data. Also, some physical reasonableness of the observations was expected. For example, when a station in a usually well-stratified region was found to have a constant density over the full depth, it was removed. To get a regular spatial representation of a density field, from which the baroclinic pressure gradient is calculated, statistical interpolation was used. It is a minimum variance estimation procedure which attempts to minimize the expected analysis error variance (see appendix F for details). In the literature, this interpolation is often referred to as optimal interpolation. It is also known by meteorologists and oceanographers as objective analysis and commonly called Gauss-Markov smoothing in oceanographic literature. 4.5.1 Density observations Only the density data obtained from cruises which covered most of the study region were chosen for this research work. Cruises done under the Special Survey in 1962, NCODE in 1984-85, and PERD in 1991-92 were considered dense enough to be used here (see Appendix B). It should be noted that the 1992 cruise was excluded from this study. As those data are now available, they could be part of a follow-up to this thesis. Chapter 4. Modelling the Dixon Entrance System z-level interval (m) 10 25 50 100 81 Position in water column (m) 0-50 50-250 250-300 300-500 Table 4.3: Interval between two consecutive z-levels. The interval depends on position in water column, the water surface being at 0 m. Cruises are divided by season. Three cruises were found for the summer time (Fig. 4.10), two for the fall (Fig. 4.11), and only one in spring and winter time (Fig. 4.12). 4.5.2 Interpolation of density data Before interpolating density data, station locations, given in earth's coordinates, had to be transformed to Cartesian coordinates. An authalic sphere, i.e. spherical model of the earth that has the same surface area as that of the actual ellipsoid, with radius of 6,371,004 m was used (Pearson II, 1990). The authalic sphere can be viewed as the intermediate step in the transformation from the spheroid to the mapping surface. In order to preserve the volumes of water, I chose to use an equal area projection, the Albers map projection with two standard parallels (Pearson II, 1990). One standard meridian and two parallels needed to be specified. These were respectively chosen as 131.75°W and, 53.5°N and 55°N. The next step was to combine any density data, which are less than 1 km away from each other, into a super-observation. This procedure facilitates the matrix inversion in the interpolation process (see Appendix F). For simplicity, the tridimensional interpolation was decomposed into 2D horizontal interpolation, using Assumption F . 6 (see Appendix F), which was executed on each of the chosen z-levels (Table 4.3). Note that a high resolution characterizes the top 50 m of water, where the pycnocline is expected to Chapter 4. Modelling the Dixon Entrance System 133W <3<W 'MW 82 I32W row 132W 133W 132W I31W 131W 131W «<!*» l » W UOW 55M 54» S31* 'MW 134W 134W |33 131 w W 133 133W W t 3 o2 2 W w 131W I 3 t w '30W 130W WOW Figure 4.10: Location (•) of C T D stations for cruises carried on in July 1-6 1985 (top), June 24-July 7 1991 (middle), and July 5-August 18 1991 (bottom). Chapter 4. Modelling the Dixon 134 W Entrance 133 W 83 System 132 W 131 W 55 rW «5N 54 B4N 131 W 134W 133 W 130 W 132 W 53 N 134W 133 W 132W 131 W 130 W Figure 4.11: Location (•) of CTD stations for cruises carried on in September 19-October 10 1962 (top) and October 14-27 1984 (bottom). Chapter 4. Modelling the Dixon 134W 134W 134W Entrance 133W 133W 133W 84 System 132 W 132W 132 W 131W 131 W 131W 130 W 130W 130 W 55 H] 54 Ki S3H 134W 133W 132W 131 W 130W Figure 4.12: Location (•) of C T D stations for cruises carried on in April 11-22 1984 (top) and December 2-18 1991 (bottom). Chapter 4. Modelling the Dixon Entrance System 85 exist. The resolution then gradually decreases to reach a minimum near the bottom. If a density value is not available at a specific z-level, a vertical linear interpolation between existing data was done. Density data were then statistically interpolated, on each z-level, from observation sites to node positions of a regular 5x5 km grid. The statistical interpolation assumes that the large-scale structure of the variable to be estimated is known (Assumption F.l). This means that we must provide the mean densityfield.As this meanfieldis unknown here, it has to be estimated from the available density data. With strong horizontal density gradients existing in the region of Dixon Entrance, the mean must be spatially variable. However, the statistical interpolation algorithm used in the thesis (called "0AX5", from Bedford Institute of Oceanography) assumes a spatially constant mean density. The algorithm was then applied locally, rather than over the whole domain, so that the impact of this constant mean is minimized. As the algorithm is applied locally, the number of nearest neighbours, representing the number of stations considered in the analysis of a specific horizontal point, has to be given. This number was set to either 5, 10, or 15 depending on the number of available stations at each z-level. 15 stations were usually chosen for the first 100 m, 10 for the next 300 m, and 5 for deeper depths. Additional information is needed before performing the interpolation. First, the variance of observation errors is needed (see Assumption F.5). As observation errors include both instrument error and error of representiveness (Assumption F.2), their variance can be considered as an adjustable parameter in the interpolation procedure. A nondimensionalized value of 0.3 was chosen. This means the ratio of the noise variance to the signal variance is 0.3. Note that a lower number would introduce physical phenomena (local density minimum or maximum) where no observations are present to confirm their presence while a higher number would smooth out the observed phenomena of interest. 86 Chapter 4. Modelling the Dixon Entrance System Secondly, the horizontal covariance function of the density variations (mesoscale component) has to be specified. Considering how few density data are available for this region, choosing a good covariance model is not an easy task. The algorithm "OAX5" automatically sets the covariance function as Cg= (l + R + ^Je- «, (4.10) R ki where Rki is defined as ^ - V—zi—+—Zj—J ' ( 4 1 1 ) L and L are the horizontal scales of correlation. These scales could also be varied with x y horizontal position but for simplicity are considered as constant over the domain. Note that Assumptions F.6 and F.7 were both applied. As the water circulation of the Dixon Entrance System is not mainly oriented in one direction, the covariance function can be further simplified by an additional assumption : A s s u m p t i o n 4.1 The horizontal covariances of the mesoscale component of density are isotropic. With Assumption 4.1, L = L = Lj and Lj is set to 15 km. (4.10) now becomes a x y two-dimensional (x,y) homogeneous isotropic covariance model which can be written as c < * - = ( i + £ + ^ ) < 4 i 2 » with Tki being the absolute horizontal distance between two points (Fig. 4.13). Examples of densityfieldsobtained with statistical interpolation at 10-m, 40-m, and 100-m depth are shown in Figs 4.14, 4.15, and 4.16. It should be pointed out that density values tend towards observation values in data-rich regions and towards the large-scale densityfieldin data voids. As the mean density is not known but estimated from limited observations, confidence in results is low in data-sparse regions. Chapter 4. Modelling the Dixon Entrance System 0 30 60 87 90 120 150 Distance Between Two Points (km) Figure 4.13: Horizontal homogeneous isotropic covariance model for density variations used in the statistical interpolation process. 4.5.3 Calculation of baroclinic pressure gradients The horizontal density gradients werefirstcalculated on each z-level using a central difference scheme and the statistically interpolated density values. The resulting gradients were then interpolated from the regular to a triangular horizontal grid using bilinear interpolation. Note that density values could have been interpolated directly from observation stations to the nodes of the triangular grid. Afiniteelement scheme would then have been used to calculate the gradients. For simplicity, the previous method was chosen. The baroclinic pressure gradient (Rrea) was calculated on each z-level using (D.13) : Rres = The integration was done following the - - f°V dz. hP Trapezoidal rule. (4.13) The next step, interpolation of Chapter 4. Modelling the Dixon 134W Entrance 133W 134W 132W 133W 134W 1 3 3 W 131 W 131 W 132W 133W 134W 88 System 132W 1 3 2 W 131W 1 3 1 W 130W 130W 130W 130W Figure 4.14: Density field, at 10-m depth, obtained from observations and statistical interpolation for June 24-July 7 1991 (top) and December 2-18 1991 (bottom). The contours are in cr units (tr =/j-1000 kg m ) with intervals of 0.1 <r . -3 t t t Chapter 4. Modelling the Dixon 134W !33W 134W 3 4 133W 133W W 89 System 133W 134W 1 Entrance 132W 132W 132W 132 W 131W 131 W 131W 131 W 130W 130W 130W 130W Figure 4.15: Density field, at 40-m depth, obtained from observations and statistical interpolation for June 24-July 7 1991 (top) and December 2-18 1991 (bottom). The contours are in <r units (o- =p-lOOQ kg m~ ) with intervals of 0.05 cr . 3 t t t Chapter 4. Modelling the Dixon 134W Entrance 133W 134W 133W 134W 133W 134 W 133W 90 System 132W 131 W 132W 131 W 132W 132W 131 W 131 W 130W 130W 130W 130W Figure 4.16: Density field, at 100-m depth, obtained from observations and statistical interpolation for June 24-July 7 1991 (top) and December 2-18 1991 (bottom). The contours are in <r units (<r=/>-1000 kg m~ ) with intervals of 0.05 a . 3 t t t Chapter 4. Modelling the Dixon Entrance System 91 Figure 4.17: Depth-averaged baroclinic forcing term (i?re») for June 24-July 7 1991. Units are cm s and values were calculated on the Extended Grid. -2 the baroclinic pressure gradients from z-levels to <r-levels, was done using vertical linear interpolation. Finally, the depth-averaged baroclinic pressure gradient was obtained by —• vertically integrating i«V , : e Xe. = J\ Rresdcr. (4.14) For comparison with the depth-averaged nonlinear tidal forcing term (Fig. 4.8), examples of Rr , are displayed on Figs 4.17 and 4.18. e Chapter 4. Modelling the Dixon Entrance System 92 Figure 4.18: Depth-averaged baroclinic forcing term {Rret) for December 2-18 1991. Units are cm s~ and values were calculated on the Extended Grid. 2 Chapter 4.6 4. Modelling the Dixon Entrance System 93 Wind and River Discharge Forcings Both the wind and river discharge effects are indirectly included in the baroclinic forcing term. The measured densityfieldis the result of an adjustment to the wind stress and the input of freshwater from rivers. These two forcings can then be viewed as local correction factors, at the free surface and at the river mouths respectively. 4.6.1 Wind stress Different types of wind observations exist for the region under investigation (see Appendix B). The ocean buoys provide the best observations since they directly give oceanic winds. Data from three buoys are available for the study region. However, they exist only since the end of the 1980's-beginning 1990's. For time periods prior to the end of the 1980's, lighthouse observations are the best data available. Note however that differences exist between winds measured onshore and offshore. Land-based wind data may not accurately reflect the wind condition over the nearby coastal waters (Schwing and Blanton, 1984; Liu et al., 1984; Hsu, 1986). Even stations located on flat, exposed land cannot represent offshore conditions (Hsu, 1981). The proper correction of the wind-velocity data over land to precisely reflect the wind stress over the nearby sea requires a knowledge of wind shear and thermal stability at both locations (Schwing and Blanton, 1984). Ng (1993) found however that a linear relationship between land-based and sea-based low-passfilteredwind-stress data is a good approximation. A transfer function can then be used to adjust the land-based wind speed so that the sea-based winds are estimated. The wind data from ocean buoys can be related to the wind data from lighthouses for common time periods. Nevertheless, since lighthouse wind observations were not available at time of data processing, the wind forcing was excluded from the model when no ocean buoy data were available. Chapter 4. Modelling the Dixon Entrance 94 System The wind velocity time series were first processed by J.Y. Cherniawsky (Institute of Ocean Sciences, Canada). Wind velocity measured at 3.7-m height was used and at 4.9-m height when the former was not available. Linear interpolation was used to fill gaps of less than two days. Since the physical quantity entering into the model is wind stress, and not the wind velocity, our interest is in wind stress. The stress was obtained from wind velocity using (3.14), (3.18), and (3.19). Examples of wind stress time series, from the ocean buoys present in the study region, are shown in Figs 4.19 and 4.20. Note that the chosen time periods correspond to the time of available density observations. Because the sea possesses inertia, the observed response of the sea is expected to depend on past wind conditions. The observed wind-driven surface current at a given time is due to the time-integrated effects of the wind stress over the time history of the wind. The main problem is here to determine the time lag of the sea response to wind forcing, i.e. choosing the time interval over which the wind forcing must be calculated in order to be consistent with the density forcing. As a first approach, the time interval corresponding to the time of the available density data is used; the time lag is set to zero. An additional simplification, spatial homogeneity of the wind field, is also made. Table 4.4 shows the time average wind speed and stress at each of the three ocean buoys with the spatially average values for two cruise periods. For comparison with the tidal —* and baroclinic forcing terms in the depth-averaged momentum equation, T , (Fig.4.8) r e and 4.6.2 lire, (Fig. 4.18) respectively, the wind forcing ($ «) is shown on Fig. 4.21. Pe R i v e r runoff River discharge is driving the model through open boundary conditions. At the mouth of both the Skeena and Nass rivers, open boundaries 11 and 13 respectively, the intensity and direction of the depth-averaged flow are specified. The current is forced to flow perpendicular to the boundary and toward the model domain. The current intensity Chapter 4. Modelling the Dixon Entrance 95 System 24 —# Figure 4.19: Stick vector diagrams of the hourly wind stress (H ^ ) for June 24-July 7 1991 at (a) buoy 145, (b) buoy 183, and (c) buoy 205 shown on Fig. B.l. Units are 10 m s with positive values for wind blowing towards the North-West and negative towards the South-East. m -5 2 -2 res Chapter 4. Modelling the Dixon Entrance 96 System 50 Figure 4.20: Stick vector diagrams of the hourly wind stress (i/ \& ,) for December 2-18 1991 at (a) buoy 145, (b) buoy 183, and (c) buoy 205 shown on Fig. B.l. Units are 10 m s with positive values for wind blowing towards the North-West and negative towards the South-East. m -5 2 -2 re Chapter 4. Modelling the Dixon Entrance Figure 4.21: Wind forcing term (* «) for December 2-18 1991. Units are cm s values were calculated on the Extended Grid. Pe 97 System 2 and Chapter 4. Modelling the Dixon Entrance 98 System Wind stress Wind velocity Intensity Direction Intensity (m s") (Degrees) (x IO" m s- ) 1 5 June 24-July Buoy 145 Buoy 183 Buoy 205 Average 2.63 0.55 1.36 1.16 Buoy 145 Buoy 183 Buoy 205 Average 2.75 5.84 5.77 4.21 103 194 50 94 December 32 337 33 11 2 2 (H ^ ) m res Direction (Degrees) 7, 1991 3.17 0.75 1.22 1.30 2-18, 1991 4.13 10.54 10.48 7.17 105 210 62 104 33 329 32 6 Table 4.4: Wind velocity and stress averaged over two cruise periods, June 24-July 7 and December 2-18 1991. The time average is given for the three ocean buoys shown on Fig. B.l as well as their spatial average. Wind direction is defined as the direction towards which the wind blows with angles measured clockwise from true north. was calculated from theflowrateof the corresponding river for the specific time of the year (see Fig. 4.22 and Table 4.5). The depth-averaged current is actually obtained by dividing the time-averaged discharge rate, calculated over the corresponding cruise period, by the width and the depth of the river mouth. Chapter 4. Modelling the Dixon Entrance 99 System Figure 4.22: Daily discharge rate of the Skeena (thin line) and the Nass (thick line) rivers for the year 1991. The vertical dashed lines show June 24, July 7, December 2, and December 18. Cruise time period July 1-6 1985 June 24-July 7.1991 July 5-August 18 1991 September 19-0ctober 10 1962 October 14-27 1984 April 11-22 1984 December 2-18 1991 Flow rate (m s ) Skeena Nass 2523 2895 1884 2386 1346 1151 631 607 728 476 489 386 378 459 3 _1 Table 4.5: Time-averaged discharge rate of the Skeena and Nass rivers computed over each of the cruise period. Chapter 5 Tidal Rectification and Wind-Driven Circulation This chapter describes the tidal residual currents arising from the interaction of tidal currents with topography (tidal rectification) and the wind-driven currents obtained under specific wind conditions. The importance of the vertical advection of tidal momentum is presented. The influence of the earth's rotation and the effect of density stratification on the currents are discussed. The sensitivity of model results to grid resolution, open boundary conditions, and friction parameters is analyzed. Comparisons with results from ear her studies and models are made throughout the chapter and a summary concludes. 5.1 Tidal Residual Currents Tidal residuals are defined here as the net flow obtained when averaging the currents over a tidal cycle (12.42 hours). These tidally-driven residuals are persistent, and thus important features. They are linked to the local bottom or coastal topography. Even if they are considerably weaker than storm-driven residual flows, they can contribute more significantly to the overall long-term distribution and transport of water properties than do the stronger, but intermittent and directionally variable wind-driven flows. Only the residualflowswhich form closed streamlines within the area being studied are simulated. Non-circulatory residualflowsdepend on the conditions imposed at the open boundaries, i.e. conditions in the outside ocean. In order to solve for these residuals, observations of mean sea-level slopes would be necessary. This study is then limited to the internal circulatory residuals. 100 Chapter 5. Tidal Rectification and Wind-Driven Circulation 5.1.1 101 Generation mechanisms The main generation mechanism of residual circulation, byfluctuatingcurrents (tides) over ridges, is the stretching and shrinking of the fluid columns (Robinson, 1983). As a reminder, the vorticity (£), also called relative vorticity, is defined as £ = dv/dx — du/dy and the potential vorticity as (/ + Q/H , where the Coriolis parameter / is named m ambient vorticity. Provided there is a component of flow normal to the topographic axis, as fluid columns move onto shallow waters, they are squeezed and the tendency for potential vorticity to be conserved results in anticyclonic (clockwise) vorticity (in northern hemisphere) being generated. Conversely as fluid columns move off shallow water into deeper water, they spin up cyclonically. We then get negative (clockwise) vorticity over the shallow area, and positive in the surrounding deep water. The interaction of the tidal currents with bottom friction constitutes the second main generation mechanism (Robinson, 1983). Consider first the case where the frictional variability, and the ensuing vorticity and residual circulation, is induced by topographic variations. Bottom friction tends to generate vorticity because the same bottom frictional stress will have more effect on a shallowfluidcolumn than on a deeper column. In other words, for the same bottom velocity, the bottom-friction force is the same, but its depthdistributed effect, through vertical eddy viscosity, is greater in shallow than it is in deep water. The bottom friction also varies with the spatial variability of the bottom current field; a fasterflowtends to be retarded more than a slower one. A third residualflowgeneration mechanism, usually neglected, is the spatial variation of bottom roughness (Gross and Werner, 1994). The contribution to the residual circulation from the gradient of bottom friction coefficient can however be as large as that arising from gradients of depth and bottom velocity. 102 Chapter 5. Tidal Rectification and Wind-Driven Circulation 5.1.2 The vorticity equation Examination of the vorticity is a useful way of understanding how the residual generation mechanisms operate. For simplicity, the dynamics of the tidalfloware examined in twodimensional (x,y) form by considering the depth-averaged equations. The total acceleration term, in the depth-averaged x-momentum equation (D.12), is expressed as 1 du ~di ~ HZ 1 du j rv dt ~HZ Z J-hm p du duu J-H, dt dx dvu dwu dy dz dz. Using Leibnitz' rule, the above expression can be rewritten as du du dt ~di = u dn {u)^ dr\ duu uu dH (TO)„ dr\ H\ ~dl ~ ~H^m ~dx~ H^~dx H ~d~x' {uu)_ dh dm vu dH (vu)n drj {vu)- dh H dx dy H dy H dy H dy m + + + n m hm m m m m hm m v m Now, introduce the kinematic boundary conditions (D.5 and D.8) and the depth-averaged continuity equation (D.ll) into (5.1) : du ti du = _du _du m o-x dy *> +U +7U +a . ( 5 - . 2 ) where expresses the x-momentum defect involved in transforming the 3D horizontalflowfield into two dimensions (Robinson, 1983). If there is any vertical structure in the horizontal flowfield,whether due to the bottom boundary friction, surface wind stress, density gradients or any other cause, terms [uu — (u) ] and [uv — (u)(y)] are generally non-zero. 2 However, it is only where there are horizontal changes in the structure, or horizontal 103 Chapter 5. Tidal Rectification and Wind-Driven Circulation gradients of total water depth, that a contributes to the total acceleration. From the x substitution of (5.2) into the depth-averaged momentum equation (D.12), we get — + u • V [u) + a + / x u + gV ri + at h (5.3) = 0, h n m with It should be pointed out that wind forcing as well as baroclinic pressure gradients are left out here because the present interest is in tidal residuals only. The vorticity equation is obtained by taking the curl of (5.3) : — = -u-V ( h - (f + ()(V -u)-V xa + h ^—•— h Cioo\ub\ub x V iir /V/iCioo C 100 fc m Vh\ub (5.4) H„ where the vorticity is denned as dv du dx dy' dv b du dx dy and the bottom vorticity as b Using the depth-averaged continuity equation ( D . l l ) , (5.4) can be rewritten as \ ( drj dt V 3 H dt m + u-V H h n Vh x a \ x Cioo \ 0 E I C\.QQ\ub\db ClQ \ub\Cb G where terms A to J are described in Table 5.1. H m H \ub\ (5.5) 104 Chapter 5. Tidal Rectification and Wind-Driven Circulation Term Physical process A B C D E F G H I Rate of change of vorticity Advection of vorticity Vorticity adjustment to accommodate surface height fluctuations Topographic stretching/shrinking of vorticity Generation of vorticity by lateral variation in vertical structure Dissipative part of bottom friction (assuming £ and of same sign) Generation of vorticity by gradient in the bottom friction coefficient Generation of vorticity by gradient in depth-distributed friction Generation of vorticity by bottom velocity gradient in bottom friction Table 5.1: Description of the terms of the vorticity equation (5.5). 5.1.3 Numerical results and analysis of vorticity terms The governing equations introduced in Chapter 3 (3.36-3.49), with modification (4.6), were solved on the Large Grid (Fig. 4.5) and 11 tr-levels (Fig. 4.7(a)). The wind forcing term (V& e») P and the baroclinic forcing terms (Rr ea and R ) were set to zero, allowing Tea one to concentrate on the tidally-driven residual currents. River runoff forcing was also excluded. The resulting surface elevation and depth-averaged currents are shown on Figs 5.1 and 5.2 respectively, and a zoom of the largest currents (geographical region delimited by the box in Fig. 5.2) on Fig. 5.3. Note that in interest of clarity, the computed velocities were interpolated onto a 3x3 km grid before being plotted. The interpolation was applied to all tidal residual velocityfieldspresented in this chapter. Also, all tidal residual currents will from now on be displayed for the zoom region (box shown on Fig. 5.2) only. The largest computed residual currents, reaching as high as 24 cm s , are found -1 just north of Langara Island (see Fig. 5.4 for the location of geographical regions discussed here). In the region surrounding Learmonth Bank, currents reach a maximum of 20 cm s _1 while along the Rose Spit Sill they reach a maximum of 15 cm s . Note that _1 Chapter 5. Tidal Rectification and Wind-Driven Circulation 105 Figure 5.1: Residual surface elevation generated by tidal forcing. The thick solid line represents the zero surface elevation contour, thin solid lines positive values and dotted-dashed lines negative values. Contour interval is 0.2 cm. 106 Chapter 5. Tidal Rectification and Wind-Driven Circulation Figure 5.2: Residual depth-averaged flow driven by tides. Units are cm s delimits the zoom region shown on Fig. 5.3. 1 and the box Chapter 5. Tidal Rectification and Wind-Driven Circulation 107 Figure 5.3: Residual depth-averaged flow driven by tides for the zoom region (box shown on Fig. 5.2). The units are cm s . -1 Chapter 5. Tidal Rectification and Wind-Driven Circulation 133W 132W 108 131 W Figure 5.4: The five geographical regions over which the terms of the vorticity equation were evaluated, with depth contours as shown on Fig. 2.1. Regions 1 to 5 correspond respectively to Rose Spit Sill, Learmonth Bank, Langara Island, Cape Chacon, and Cape Muzon. 109 Chapter 5. Tidal Rectification and Wind-Driven Circulation the tidal residual currents computed along this sill are of the same order of magnitude as the measured currents (6-22 cm s ) (Bowman et al., 1992) but the Rose Spit Eddy Section 2.5) is not fully formed here. However, the Gyre is fully formed _1 Learmonth Bank (see and, from drifter simulations performed at 10-m depth, characterized by a rotation period of about 3 days. It should nevertheless be pointed out that none of the numerical particles deployed in the surrounding of the gyre completed more than one loop. Most of them were instead trapped into the clockwise eddy located north of Langara Island (see Fig. 5.3). As matter of interest, the 3D particle tracking algorithm for afiniteelement grid with linearfiniteelements, written by Brian 0. Blanton from University of North Carolina in 1995, was used for all drifter simulations performed for this thesis. As seen previously, the vorticity equation can help us understand how the residual currents are generated. First, introduce the harmonic expansion in time (D.14-D.16) in equation (5.5), assume a spatially-constant bottom friction parameter (see Section 4.3), and look for the steady-state solution : o= - T 22 P • U V ^ - P + * p=-N P*O 77 _ v , m +1 £ BB . L VfcJ^m "77—l re> 7jr~TC- « e ^ , P ™ ' EE CC -4, u ™>r a m AA ClOO x — KJ JE * p=-N *o V/i7 C100 —TT— + x -77-7«re. x FF /c e\ + higher-order terms, (5.6) GG where N a. 4# ^ ( Q _ I dx~ { ( ^ W - p V p ~ (*-P^»))] m + ^- [{h )(u-pV m p - {u-p){vp))]^ + higher-order terms. Table 5.2 provides the physical description of terms A A to GG. Note that weak nonlinearity was applied, i.e. the interaction terms of the residualflowwith itself are neglected, and N is the number of tidal constituents (see Table 4.2 for a list of those used here). Chapter 5. Tidal Rectification and Wind-Driven Circulation Term AA BB CC DD EE FF GG 110 Physical process Advection of tidal vorticity by the tides Topographic stretching/shrinking of residual vorticity Topographic stretching/shrinking of tidal vorticity Generation of vorticity by lateral variation in tidal vertical structure Dissipative part of bottom friction Vorticity generated by depth gradients in depth-distributed friction Vorticity generated by bottom velocity gradient in bottom friction Table 5.2: Description of the terms of the steady-state vorticity equation (5.6). Term AA BB CC DD EE FF GG Region 1 Region 2 Region 3 Region 4 Region 5 4.254 2.913 3.898 4.050 2.361 1.535 2.027 1.476 2.140 1.961 6.214 2.022 0.728 1.573 0.553 7.442 6.148 1.012 2.001 1.463 3.938 0.571 0.925 0.376 0.398 4.343 0.355 0.217 0.529 0.285 0.122 1.479 0.226 0.061 0.081 Table 5.3: Absolute value of each term of the steady-state vorticity equation (5.6) averaged over the five regions shown on Fig. 5.4. Units are 10 s . _lo -2 Each term of the vorticity equation (5.6) was evaluated over the five geographical regions shown on Fig. 5.4, regions characterized by important tidal residual currents. The spatial average of the absolute value of each vorticity term is given in Table 5.3. In region 1, namely Rose Spit Sill, the most important generation mechanisms of residual vorticity, in order of importance, are the stretching/shrinking of thefluidcolumns by the tides, the interaction of the tidal currents with bottom friction (through topographic variations), and the advection of tidal vorticity by the tides. In region 2, around Learmonth Bank, advection is the most important mechanism with topographic stretching/shrinking of residual and tidal vorticities next. In both regions 3 and 4, close to Langara Island and Chapter 5. Tidal Rectification and Wind-Driven Circulation 111 Cape Chacon, residual vorticity is mainly generated by lateral variations in the tidal vertical structure, with advection of tidal vorticity as the second mechanism. Finally in region 5, namely Cape Muzon, the advection of tidal vorticity and the topographic stretching/shrinking of residual vorticity are, in order of importance, the main generation mechanisms. 5.2 Role of the Vertical Advection of Tidal Momentum by the Tides The numerical estimates of steady-state residual velocities presented in Ballantyne et al. (1996) for the Dixon Entrance region were obtained assuming that vertical advection was negligible. This section discusses the importance of vertical advection through momentum and vorticity analyses, as well as numerical results. 5.2.1 The momentum equation The advection term in the x-momentum equation (D.2) is expressed by duu duv dx „ du dy du duw dz dv du d~z dw (5.7) and the continuity equation (D.4) as (5.8) Introducing (5.8) into (5.7), we get that _ . _r*. Ou du Ou V • (utl) = u——h v——\w dz~~ (5.9) If we now neglect the vertical advection term (^ ), we get L _ , f*. „ Ou V.(«t/)=2» g ; + Ou Ov ^ ^ , + (5.10) Chapter 5. Tidal Rectification and Wind-Driven Circulation 112 or, introducing (5.8), _ . ->. du T V. du dw = (5.11) Neglecting the vertical advection is then equivalent to substituting wj^ by i n the momentum equation. For the region under study, the effect of this substitution on the depth-averaged tidal advection term ( T , ) is displayed on F i g . 5.5. r e 5.2.2 T h e vorticity equation If the vertical advection is neglected, (5.2) then becomes du ti du _du m dx- = +u + _du dy~ 7U ° + a (5.12) *> + b with i \ dn / . dn , . (u)„-— + (uu)„——h [UU)-L— b = m dh dn dh h [vu)„——h (VU)_L m m —— Now, equation (5.3) is replaced by du + u • Vfc(«) + a + 6 + / x i i +flrVfcr?+ = 0, at (5.13) n m where b= + (-5)_ [(iZ)_ .V fc ] fcm Hr, hm fc m and equation (5.5) by ( Jjt, A (f + C) dn + H m B - u-V H„ v Cioo\Ub\(b Vixa- h •> — Hr, E D ( ClOO\Ub\Hb ^ H m VhCioo V H h Cioo m H m H | V \u \ h b y^x_6, \ub\ I / J (5.14) Chapter 5. Tidal Rectification and Wind-Driven Circulation 134W 133W -•. ^ res - 2 132 W 131 W 132W 131 W •55 N .0.0005 . 0.001 133W Figure 5.5: momentum f (x,y) cm s and 113 132 W -64 N 131 W Difference between the depth-averaged nonlinear tidal forcing term in the equation when vertical advection is respectively neglected and included, f ,(x,y,z), (top) and f ,(x,y,z), zoom of Fig. 4.8, (bottom). Units are values were calculated on the Large Grid. re re Chapter 5. Tidal Rectification and Wind-Driven Circulation 114 Region 1 Region 2 Region 3 Region 4 Region 5 2.001 1.117 3.297 3.178 3.873 69 18 85 43 63 Term HH % Table 5.4: Absolute value of term HH of the steady-state vorticity equation (5.15) averaged over the five regions shown on Fig. 5.4. Units are 10 s . The % represents the relative importance of term HH with respect to the most important vorticity term (bold value in Table 5.3). -lo -2 where term J represents the extra generation of vorticity introduced by the absence of vertical advection. When the harmonic expansion in time is introduced, the steady-state form of (5.14) becomes, similarly to (5.6), 0= - l £ „ , fu .-V H re h , 1 ^ m (p"--p u V/jHjn „ _ -Vixa, H p n p^O p#0 BB AA ClOO "7^re« H„ Cioo DD CC VhH u H„ C o -4 1" ~ 7 J — m x — Vfc7 1 0 H„ FF "fc x GG x £ I . /c i c \ °res + higher-order terms, ( 5 . 1 5 ) HH with 1 Tjf m N £p=-N *5 p?0 ft-p ' ^hh ) m + higher-order terms. The only difference between the vorticity equations (5.6) and (5.15) is the extra term HH present in the latter. Excluding the vertical advection of tidal momentum by the tides has then the effect of generating extra vorticity and associated residual currents through term HH. Term HH was evaluated for the same geographical regions over which the other terms of the vorticity equation were calculated (Fig. 5.4). The resulting values are given in Table 5.4. Assuming that vertical advection is negligible then corresponds to adding very Chapter 5. Tidal Rectification and Wind-Driven Circulation Figure 5.6: Difference between the residual depth-averaged flow driven by tides when neglecting the vertical advection and the residual depth-averaged flow driven by tides with vertical advection included. The units are cm s . -1 large residual vorticity around Learmonth Bank, large residual vorticity in the neighbourhood of Capes Chacon and Muzon, considerable residual vorticity north of Langara Island, and small residual vorticity along the Rose Spit Sill. 5.2.3 N u m e r i c a l results The impact, on the residual depth-averaged currents driven by tides, of neglecting the vertical advection is shown on Fig. 5.6. As expected from the vorticity analysis, the 115 Chapter 5. Tidal Rectification and Wind-Driven Circulation 116 most noticeable difference is around Learmonth Bank. Note that the inclusion of vertical advection does not really alter the circulation pattern but decreases the intensity of the flow. Without the vertical advection, residual velocities are higher by as much as 48 cm s _1 near Learmonth Bank. Ballantyne et al. (1996) solved the problem of modelling excessively large residual currents, compared to observations, by adding extra smoothing to the nonlinear tidal forcing terms of the momentum equation (T Pe4 and T ). This smoothing had the effect rea of decreasing the tidal forcing terms, therefore removing the additional tidal advection due to the exclusion of vertical advection, but also the effect of spreading out the forcing terms and therefore the interesting circulation patterns. The anticyclonic eddy obtained by Ballantyne et al. (1996) around Learmonth Bank is spatially larger than the eddy obtained in this thesis. I also studied the impact, on the tidally-rectified currents, of including or not the vertical advection for other geographical regions. For example, the flow around Cape St.James (Fig. 1.1), which is the most dynamic and complexflowof the north BC coast (Crawford et al., 1995), was studied. In this region, the maximum speed obtained when excluding the vertical advection of tidal momentum is 44 cm s . Including the vertical -1 advection decreases the maximum speed to 30 cm s , a value closer to the Crawford et _1 al. (1995) estimate, from observations, of 24 cm s . Note that the barotropic version of _1 the model used here gave a maximum value of 69 cm s _1 (Foreman et al., 1993) while Ballantyne et al. (1996) calculated a maximum speed of 59 cm s . _1 In conclusion, the inclusion of vertical advection of tidal momentum by tides is necessary because it is important, as shown by the fact that it leads to speeds much closer to those observed. Chapter 5. Tidal Rectification and Wind-Driven Circulation 5.3 117 Wind-Driven Currents We have up to this point emphasized the importance of tidal rectification in driving the meanflowin the Dixon Entrance System. There is no doubt that the winds have some influence upon the currents. Steady-state wind-driven currents were computed the same way as the tidal-rectified currents. The only difference lies in the forcing terms; the nonlinear tidal forcing terms Q s, re ^res, and T are now set to zero and the wind res forcing term \& , turned on. re 5.3.1 Idealized wind Before discussing the influence of the real wind, we study the currents generated by idealized winds. These simulations provide some insight into the dynamics of the system before applying real conditions. The results are compared with the depth-averaged wind currents obtained by Hannah (1992) in his barotropic model. Four spatially uniform steady winds were used for this study: 1. Wind blowing from the South-East at 10 m s 2. Wind blowing from the East at 10 m s _1 _1 3. Wind blowing from the North-West at 10 m s 4. Wind blowing from the West at 10 m s -1 _1 The first and second wind directions correspond to winter conditions in Hecate Strait and Dixon Entrance respectively while the third and fourth wind directions correspond to summer conditions in HS and DE respectively (see Section 2.3). For comparison purposes, the same wind speed was used for the four wind directions. Figures 5.7 and 5.8 display the current vectorfield,at 10-m depth, driven by a wind blowing from the SE and a wind blowing from the East respectively. Notice that the 118 Chapter 5. Tidal Rectification and Wind-Driven Circulation 134W 133W 132 W 131 W 130 W * t T . . f r- r-, , , 1 1 1 , . . . . , p * * P* f r p p p r p P p . ***/ *T*pppp*P*. P ****** pp pp * p r . t 55 P * * * P * * P P P P * P P P **?**pppp*tpp * * * * * * * * * * * ****** f r * * - Tg p .f, r * * * *%f ***T*P*.>i * * P* * * \ * t * * * * * P r i t t r t T . * * * * T * P * * * t T T t 1 r . . * * * * * * * 7 < * 7 > * * A * » * P . . r * * * * * * P * * * * * * * * * * * p *t*** * p p * * * * ****** * * * * * * * * * * * p******* * * * * ***7>7>P . P PP r * * r r ll^'^'i *******»yi*Ji»*7'*?** P P * *7<**P P * * ****SiS'S'S<S'Sl}'* ***??****?*' p * * * * * * * * * * * * * * * i * * * * * * * * * * * * P P P p *** * * ?j**j<j;siMS'?'** ************* p *** * ? 7> * /* PJLTCi. *?* * ' t * * r t t r 1 P t P * t ? 7> f ? -ty 54 N- 134 W ZCtS. » f » » » » » » t t t f » t « 133W 132 W 131 W 130W Figure 5.7: Steady-state currents, at 10-m depth , driven by a SE wind of 10 m s . The 1 - l units are cm s Chapter 5. Tidal Rectification and Wind-Driven Circulation 134W 133W i — 134W 1 133W 131 W 132W -— ^ r* 132W =- 1 131 W 119 130W 1 130W Figure 5.8: Steady-state currents, at 10-m depth, driven by a easterly wind of 10 m s The units are cm s . _1 Chapter 5. Tidal Rectification and Wind-Driven Circulation 120 wind-drivenflowpattern obtained with a wind blowing from the NW and from the West can simply be found by reversing the vectors in Figs 5.7 and 5.8 respectively. This is actually the case because the model dynamics are linear here. The model used in this thesis is weakly non-linear (see Section 3.5) and then becomes linear whenever the nonlinear tidal forcing is absent. Note also that all wind-driven currents presented in this chapter were interpolated on either a 5x5 or 6x6 km grid before being plotted. To compare model results with Hannah's results (1992) (Fig. 5.9), the depth-averaged flow for a SE wind is shown in Fig. 5.10. Large velocities are present in the deep channel of Hecate Strait (east side), in the channel connecting HS with Chatham Sound and also in the sound. This confirms the presence of the Chatham Sound Diversion found by Hannah (1992), i.e. the diversion, under SE winds, through Chatham Sound of a significant fraction of the water transported from Hecate Strait to Dixon Entrance. Note also that the wind-driven circulation simulated here, and the current vector field computed by Hannah (1992) as well, shows no evidence of any eddy in Dixon Entrance. The major difference between Figs 5.9 and 5.10 lies in the currents on the west coast of the Queen Charlotte Islands. Note however that Hannah's results did not reach a steady state in the deep ocean. Another remarkable difference is the opposite direction of the circulation in the northern half of Dixon Entrance. 5.3.2 Observed wind Figure 5.11 shows the wind-driven currents, at 10-m depth, forced by the wind forcing term shown on Fig. 4.21. This velocityfieldcorresponds to winter conditions, specifically to the average wind condition of December 2-18 1991. The maximum speed, 16 cm s , _1 is found in Chatham Sound. From Figs 5.2 and 5.11, it can be concluded that the surface currents driven by a ~ 2-week average winter wind forcing are more important, in some areas of the Dixon Entrance System, than the tidal residual currents. It should be kept Chapter 5. Tidal Rectification and Wind-Driven Circulation 121 Figure 5.9: Depth-averaged wind currents for a SE wind obtained by Hannah (1992). 122 Chapter 5. Tidal Rectification and Wind-Driven Circulation 134 W i 134W 133 W . , 133W 132 W r*- 132W 131 W 1 131 W Figure 5.10: Depth-averaged wind currents for a SE wind of 10 m s cm s . _1 130 W 1 130W The units axe Chapter 5. Tidal Rectification and Wind-Driven Circulation 123 Figure 5.11: Current vector field, at 10-m depth, driven by the average wind observed in the December 2-18 1991 time interval. The units are cm s" . 1 124 Chapter 5. Tidal Rectification and Wind-Driven Circulation in mind that the impact of the wind on the water circulation decreases with water depth. The same analysis was done for the wind forcing term corresponding to the average wind condition calculated over the June 24-July 7 1991 time period (see Table 4.4). The maximum speed then obtained at 10-m depth is 1 cm s _1 in Chatham Sound. The wind- driven currents arising from a ~ 2-week average summer wind condition are therefore small compared to the tidally-rectified currents over most of the study area. 5.4 Influence of the Earth's Rotation To test the sensitivity of the model results to the Coriolis force, the numerical model was rerun but without rotation (/ = 0). In this way, model results can be compared with the non-rotating Hecate Model (Bell and Boston, 1962; 1963) and Hannah's results (Hannah, 1992). The depth-averaged tidally-driven residualflowfor a non-rotating earth is shown on Fig. 5.12. The first noticeable difference with Fig. 5.3 is an increase in speed. The velocities north of Langara Island now go up to 60 cm s , around Learmonth Bank up _1 to 45 cm s , and along Rose Spit Sill up to 25 cm s . -1 _1 The major difference in the circulation pattern is the disappearance of the eddy around Learmonth Bank. Note that the bulk pattern of the tidal residuals along the Rose Spit Sill remains the same, being consistent with results from the Hecate Model. Notice that the high velocities obtained with the Hecate Model, as high as 50 cm s _1 north of Rose Spit, can be partially explained by the absence of rotation in the model. As suggested by Crawford and Thomson (1991), those strong residual currents may have also been partially generated by the magnification of the frictional forces present in a hydraulic model. The depth-averaged wind-drivenflow,by a SE wind of 10 m s , for a non-rotating _1 earth is shown on Fig. 5.13. As for the tidal residuals, the non-rotating model gives 125 Chapter 5. Tidal Rectification and Wind-Driven Circulation 54 N 54 134 W 133W 132 W 131 W N Figure 5.12: Depth-averaged tidal residual velocities on a non-rotating earth. The units are cm s . _1 Chapter 5. Tidal Rectification and Wind-Driven Circulation 134W 133W 132 W 126 130 W 131 W 54 N- 53 N 134W 133W 132 W Figure 5.13: Depth-averaged velocity field for a SE wind of 10 m s earth. The units are cm s . -1 130W 131 W 1 and a non-rotating Chapter 5. Tidal Rectification and Wind-Driven Circulation 127 higher velocities almost everywhere in the domain (see Fig. 5.10 for comparison). The quite intense current on the north of Langara Island has however disappeared and the current direction of northern Dixon Entrance has reversed. Hannah (1992) also studied the influence of the earth's rotation on the wind-driven circulation (see Fig. 5.14). He obtained a circulation pattern similar to the one presented here. 5.5 Sensitivity of the Model Results In order to have confidence in the model results, one has first to make certain that the open boundaries do not greatly affect the numerical solution in the region of interest. Secondly, the grid resolution has to be fine enough to avoid numerical errors. There Eire however two limiting factors on how fine a grid can be: computer time increases with refinement and we are limited by the known resolution of the bathymetric map. Finally, as the values of the friction parameters are known only within a certain range, the numerical solution should not be too sensitive to the specific choice of the values inside that range. 5.5.1 Open boundary conditions The best way to avoid specifying any open boundary conditions would be to use the Large Grid (Fig. 4.5) for all simulations. The problem is that density observations are only available within a restricted region (see Figs 4.10-4.12). In the perspective of simulating the complete circulation driven by the four forcings, we then have to make sure the open boundaries of the grid to be used, the Extended Grid (Fig. 4.4), will not affect the solution. To study the effect of open boundaries on the tidal residual currents, the model was run on both the Large and Extended grids. There are some important differences Chapter 5. Tidal Rectification and Wind-Driven Circulation Figure 5.14: Depth-averaged wind-driven flow for a SE wind and a non-rotating obtained by Hannah (1992). 128 Chapter 5. Tidal Rectification and Wind-Driven Circulation 129 in the solution along some open boundaries of the Extended Grid but the effect of the boundaries on the interior solution is negligible. Note that most of the tidally-driven flow of the Dixon Entrance System is generated inside this geographical region itself (see Fig. 5.2). There is no importantflowinto or out of the model domain. Therefore, setting the open boundaries closer to Dixon Entrance does not alter the internal residual circulation driven by the tides. On another hand, the open boundaries have a non-negligible effect on the wind-driven circulation of the region under study. The difference in model results, obtained using the Large and Extended grids, is shown on Fig. 5.15 (see Fig. 5.11 for comparison). The current vectorfieldobtained when running the model on the Extended Grid is missing the circulation pattern due to the water exchange of the Dixon Entrance System with the outside ocean. One way to overcome the problem is to specify the inflow/outflow going through the open boundaries 1 (Revillagigedo Channel), 2 (Clarence Strait), and 7 (Hecate Strait) (see Fig. 4.4). The depth-averaged velocities along these open boundaries can be extracted from the solution calculated on the Large Grid. If we rerun the model on the Extended Grid with open boundary conditions extracted from the Large Grid solution, we now get a negligible difference between the results obtained on the two grids. The only noticeable differences are located close to the open boundaries and don't affect the solution of interest. 5.5.2 G r i d resolution There are two ways to improve the grid resolution, that is to say horizontally and vertically. The model sensitivity to both horizontal and vertical grid resolutions is analyzed here. To study the impact of the horizontal resolution, the tidally-rectified and winddriven currents were both computed on the Coarse (Fig. 4.2) and Fine (Fig. 4.3) grids using 11-sigma levels. To study the impact of the vertical resolution, the currents were Chapter 5. Tidal Rectification and Wind-Driven Circulation 130 Figure 5.15: Difference between the 10-m depth current vector field, driven by the average wind measured over December 2-18 1991, computed on the Large Grid and on the Extended Grid. Chapter 5. Tidal Rectification and Wind-Driven Circulation 131 54 Figure 5.16: Difference between the depth-averaged tidal residual currents computed on the Fine and the Coarse grids. Units are cm s . -1 computed on the Extended Grid and both 11- and 21-sigma levels (Fig. 4.7). Figure 5.16 displays the effect of the horizontal grid resolution on the depth-averaged tidal residual currents (see Fig. 5.3 for comparison). The solution computed on the Coarse Grid is definitively missing some important features related to the bathymetry. Therefore, the currents should be modelled on the Fine Grid or a grid with similar resolution, such as the Extended Grid. Note however that the refinement process is not completed here. The grid should actually be refined until no extra features appear. As further grid refinement is beyond this thesis, it should then be part of a follow-up. Note Chapter 5. Tidal Rectification and Wind-Driven Circulation 132 Figure 5.17: Tidal residual velocities in Dixon Entrance obtained by Ballantyne et al. (1996). Single shafts in multi-shafted vectors represent 5 cm s . -1 also that numerical results presented in Ballantyne et al. (1996) were obtained using what we call the Coarse Grid. It is then not surprising to see a different circulation pattern in the neighbourhood of Langara Island (see Figs 5.3 and 5.17). The horizontal resolution of the grid has much lesser effect on the wind-driven currents. It should be pointed out here that a fully-formed Rose Spit Eddy is simulated on the Coarse Grid. From drifter simulations performed on my model results, a rotation period of about 100 days was found. This is similar to what Ballantyne et al. (1996) expected from the current vector field they computed but much longer than what is expected from Chapter 5. Tidal Rectification and Wind-Driven Circulation 133 Figure 5.18: Difference between the depth-averaged tidal residual currents computed with 11 and 21 vertical sigma-levels. Units are cm s . _1 the observations (20-30 days (Crawford and Greisman, 1987) and 5-40 days (Crawford, 1995)). As this fully-formed gyre is not observed in the velocity field obtained on the Fine Grid, or on the Large Grid (see Fig 5.3), it then appears that its presence is only an artifact of the low horizontal grid resolution. Indeed, tidal residual circulation can not generate, by itself, a fully-formed Rose Spit Eddy. This result supports the conclusion made by Bowman et al. (1992) that the eddy generated by the Hecate Model (see Fig. 2.13) could have actually been generated by phase errors and then simply be an artifact. Figure 5.18 displays the effect of the vertical grid resolution on the depth-averaged Chapter 5. Tidal Rectification and Wind-Driven 134 Circulation tidally-rectified currents. There are some differences in the current intensity (maximum speed difference of 7 cm s _1 around Learmonth Bank) but not really in the circulation pattern. Considering the numerical error arising from other sources, such as we will see for the friction parameterization, the vertical grid composed of 11 elements seems to be fine enough for our purpose. Note that the effect obtained on the wind-driven currents is even smaller than the one for the tidal residuals. 5.5.3 Friction parameters A sensitivity test to friction parameters was investigated in order to evaluate the impact of possible errors in the friction parameterization of the model. All results were computed on the Extended Grid and 11 vertical sigma levels. First, the bottom friction coefficient (Cioo) was decreased from 0.01 to the value used by Naimie et al. (1994), i.e. 0.005 (see Fig. 5.19). As it might be expected, decreasing the bottom friction enhanced the currents. The maximum speed increase, 12 cm s , is _1 located in the neighbourhood of Learmonth Bank. Doubling the coefficient value has a similar effect but in the opposite sense, i.e. retarding the currents. It should be pointed out here that the influence of the bottom friction coefficient value on the tidal residual currents is more important than the influence of the vertical resolution obtained in the previous subsection. Wind-driven currents are much less sensitive to the value of the bottom friction coefficient than tidal residuals. Based on its value range (see Section 3.1), the coefficient for bottom and surface boundary layers (Ci) was increased from the value 0.3 to 0.5. This increase has the effect of expanding both boundary layers. If the boundary layers are thicker, the friction is felt farther away in the water column and then the currents are expected to be retarded. See Fig. 5.20 for such impact on the tidally-driven currents. If C\ is decreased to 0.2, the currents are enhanced by a similar amount. Note that in this case, wind-driven currents Chapter 5. Tidal Rectification and Wind-Driven Circulation 135 Figure 5.19: Difference between the depth-averaged tidal residual currents calculated with Cioo=0.005 and L7 =0.01. Units are cm s . -1 100 Chapter 5. Tidal Rectification and Wind-Driven Circulation 136 Figure 5.20: Difference between the depth-averaged tidal residual currents calculated with Ci=0.5 and Ci=0.3. Units are cm s . _1 Chapter 5. Tidal Rectification and Wind-Driven Circulation 137 are influenced as much as are tidal residuals. 5.6 Stratification Effect Up to this point, we have neglected the effect of stratification on residual currents. As seen in Chapter 3, when a water column is stably stratified, the vertical mixing is reduced and so is the vertical eddy viscosity, implying less friction in the water column. The modification brought by the inclusion of stratification to the governing equations is first described. The impact on both tidal residual and wind-driven currents is then presented. 5.6.1 Modification to the governing equations The effect of stratification on vertical mixing is mathematically represented, as seen in Chapter 3 (equation (3.34)), by A = F(Ri)A° , v v (5.16) where the superscript "0" denotes the case of neutral stability. F(Ri) is defined, from (3.32), (3.33), and (3.35), as F(Ri) with 1 (1 + XiRi)*'' g Ri (5.17) Apo p» dz or, with the introduction of the harmonic expansion in time (D.15-D.16) a PO 8 Ri = -r p» dz Note that X\ and X are two adjustable parameters still to be determined. 2 (5.18) 138 Chapter 5. Tidal Rectification and Wind-Driven Circulation A range of values for the stratification parameters X\ and X can be found in the 2 literature. For example, Munk and Anderson (1948) use X\ = 10 and X = 0.5, Bowden 2 and Hamilton (1975) use Xi = 7 and X — 0.25, and Pacanowski and Philander (1981) 2 use X\ — 5 and X — 2. Following the procedure introduced by Naimie et al. (1994), we 2 use the Munk and Anderson definition. Equation (5.17) then becomes F(Ri) = . ' (1 + lOiK)* (5.19) ' r V V F(Ri), and consequently Ri, must satisfy some limiting conditions. As we want to model the reduction of turbulence by the stratification, A v must be bounded by the values 0 and A°. Using (5.16) and (5.19), we have the hmiting conditions A v ->• 0 for F(Ri) ->• 0 and therefore Ri ->• oo A „ -> A° for F(Ri) 1 and therefore Ri ->• 0, which sets Ri as a positive number. Note that if Ri < 0, density variations would enhance the turbulence. This case is met when density inversions are present in the water column. One has to be careful because density inversions can be created by the statistical interpolation process. Regions of data voids and large bathymetric gradients are subject to this artifact. In order to avoid a non-physical increase of turbulence, a lower limit was put on Ri : if Ri < 0 then Ri = 0. (5.20) Only two terms of the governing equations (3.36-3.49) involve the eddy viscosity coefficient. Terms (3.45) and (3.48) are then the only ones that need to be modified in order to introduce the stratification effect into the model. They respectively become, using modification (4.6), Al = F(Ri) [/so(C o)»/6(*)7 + *o (H \%es\)' 10 m l.(z) + A*] , (5.21) 139 Chapter 5. Tidal Rectification and Wind-Driven Circulation and T rea Jn^Ah ^(Cioo) , , x JL -A K = F(Ri) . 2 sFlit du- f (5.22) with F(Ri) denned by (5.19), (5.18), and (5.20). 5.6.2 Numerical results Contrary to what Naimie et al. (1994) found, the influence of density stratification on the tidally-rectified currents is negligibly small here. This result comes from the fact that the stratification is important in the surface layer while the vertical eddy viscosity coefficient is important in the bottom boundary layer only. Actually, outside the bottom boundary layer, the neutral viscosity coefficient is equal to the background eddy viscosity (A° — A* ). A is small; decreasing its value does not make a big difference on the final 8 b9 solution. The parameterization for the friction in a water column used by Naimie et al. (1994) differs from the one used in this thesis. They define the vertical eddy viscosity as (5.23) A = F{Ri)N \^\ + A J>, 2 v b Q where A^o is set to 0.2 s, as opposed to our definition A = F(Ri) [ACo(<7ioo)*f&(z)|u | + A ] . h a v D v (5.24) Note that formulation (5.23) is only good for deep waters, where the bottom boundary layer thickness is rotationally limited (Davies and Jones, 1995). Also, the magnitude of the eddy viscosity used in (5.23) is determined by the depth-mean current. A physically more appealing velocity scale is the bed friction velocity u\ (Davies, 1993) which is used in our formulation (|-u^| = (Cioo)^|""£>|)- I Naimie et al. (1994), the neutral vertical n viscosity is constant over the entire water column and A then varies vertically with v Chapter 5. Tidal Rectification and Wind-Driven Circulation 140 density stratification only. The stratification is therefore expected to have an important impact on the currents, as obtained. In this thesis, A varies, in the neutral case, over the v water column (increasing linearly from the bottom up to a maximum and then decreasing to a minimum outside the bottom boundary layer). As for the tidal residual currents, the influence of the stratification on the winter winddriven currents is negligibly small. The eddy viscosity coefficient has now a high value in both bottom and surface boundary layers. However, as seen previously, the weak density stratification in the bottom boundary layer leads to a small Ri and therefore near-neutral stability (F(Ri) close to 1). In the surface boundary layer, the density stratification is usually large but so is the mechanical mixing, due to winter wind stress. Indeed, the stratification effect is cancelled by the decrease of wind-driven current speeds with depth. Note that during summer time the intensified density stratification is expected to decrease the wind driven currents but as these are already small, the effect is negligible. 5.7 Summary The model simulations, with momentum and vorticity analyses, showed that : • Important tidal rectification is present 1. north of Langara Island, due mainly to the generation of vorticity by lateral variation in tidal vertical structure, 2. around Learmonth Bank, due mainly to advection of tidal vorticity by the tides, 3. along Rose Spit Sill, due mainly to topographic stretching/shrinking of tidal vorticity, Chapter 5. Tidal Rectification and Wind-Driven Circulation 141 4. in the region of Cape Chacon, due mainly to the generation of vorticity by lateral variation in tidal vertical structure, 5. and in the region of Cape Muzon, due mainly to advection of tidal vorticity by the tides. • The vertical advection of tidal momentum by the tides should not be neglected. The inclusion of this term decreases, mostly in Learmonth Bank region, the tidal residual speeds to more realistic values. • The wind-driven currents, generated from the observed wind forcing averaged over about a two weeks period, are small compared to tidally-rectified currents during summer time b u t not during winter time. • The rotation of the earth 1. slows down the residual currents, driven by both tides and wind, 2. is essential for the formation of an anti-cyclonic gyre, driven by tides, around Learmonth Bank, 3. is not essential for occurrence of tidal rectification along the Rose Spit Sill, 4. and reverses the direction of the wind-driven circulation in northern Dixon Entrance. • The model results are sensitive to open boundaries when driven by the wind but not when driven by tides. The inflow/outflow needs to be specified along open boundaries of the Extended Grid to correctly model the wind-driven circulation. • The tidal residual circulation is really sensitive to the horizontal resolution of the grid (bathymetric resolution). The effect on the wind-driven circulation is negligible, as is the sensitivity of both circulations to the vertical resolution. A grid Chapter 5. Tidal Rectification and Wind-Driven made of the Circulation 142 Extended Grid in the horizontal and 11-sigma levels in the vertical gives satisfactory solutions. • T h e model results are sensitive, to a certain degree, to the value of the bottom friction coefficient when driven by tides but not when driven by the wind. The sensitivity of the model to the coefficient for bottom and surface boundary layers is less. • T h e stratification effect on both tidally- and wind-driven currents is negligible. T h e depth-averaged tidal residual circulation is found on Figs 5.2 and 5.3 and the winter wind-driven circulation, at 10-m depth, on F i g . 5.11. Chapter 6 Currents Arising from River Runoff and Baroclinic Effects The water circulation driven by river runoff input and density differences is described in this chapter. Water currents resulting from river runoff effects, through open boundary conditions, are first discussed. Next, buoyancy currents arising from baroclinic pressure gradients characterizing each season of the year are presented. The sensitivity of model results to open boundaries, grid resolution, friction parameters, as well as the influence of the earth's rotation and of the stratification on the eddy viscosity are analyzed. The chapter ends with a short summary. 6.1 Currents Driven by River Runoff It should be recalled that most of the effect of freshwater runoff is indirectly included in the baroclinic forcing term through the density distribution. The currents discussed here are driven by the river runoff momentum forcing which may be seen as a local correction factor applied at river mouths. To study the impact of this forcing on the circulation field, the numerical model was forced with river discharge, through open boundary conditions (see Section 4.6), only. The forcing terms Qre*, T ,, T , re Pe4 * P , , R\ >, Pe e and R r e , were all set to zero. The driving force used is the time-averaged river discharge measured over the July 1-6 1985 time period. It corresponds to the highest average riverflowrateamong the seven cruise periods considered in this thesis (see Table 4.5). Thus, the effect of river runoff for the other cruise periods would be smaller than the one discussed here. The computed currents, at 10-m depth, reach a maximum value of 5 cm s 143 _1 at the Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 144 Skeena River mouth (see Fig. 2.1 for the geographical location of rivers). As this region is characterized by shallow waters and high tidal turbulence (see Section 2.5), energy is rapidly (in space) removed and current speeds decrease away from the river mouth. The Nass River region is characterized by deeper waters, depths over 300 m, and therefore smaller currents. Note that speeds higher than 0.3 cm s _1 are only found close to the Skeena and Nass rivers. In most of the Dixon Entrance System, the currents hardly reach 0.1 cm s , whereas the speed of tidal rectified currents obtained in the previous chapter _1 is over 1 cm s . Consequently, the currents driven by river runoff are negligibly small _1 compared to tidal rectification, with the areas immediately near the Skeena and Nass river mouths as exceptions. 6.2 Currents Driven by Baroclinic Pressure Gradients —• —* The driving forces, represented by Rres and its depth-averaged value R r e a , were calculated from available density observations (see Section 4.5). It is the only driving mechanism —* —¥ —* introduced into the model here. The nonlinear tidal forcing terms Q ,, T , , and T , , —• re re re the wind forcing term ^ V e * , and the river runoff forcing were all excluded. Seven different baroclinic forcings, resulting from seven cruises (see Figs 4.10 to 4.12 for the location of the stations), were used. As Rret contains horizontal derivatives of density, information is needed over a wide area. Density coverage limits the extent of the horizontal grid that may be used: the model domain must be smaller than the region covered by the observations. Therefore, it is not possible to use, as we did for the tidal residual and wind-driven currents, the Large Grid (Fig. 4.5) to either avoid dealing with the open boundaries or help setting the open boundary conditions on a smaller grid. Nevertheless, the open boundaries should be kept as far as possible from areas of interest. Consequently, the grid used to Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 145 compute the buoyancy currents is the Extended Grid (Fig. 4.4) with 11 vertical sigma levels (Fig. 4.7(a)). Note that all model results arefirstinterpolated onto a regular grid (4x4 km) before being plotted. 6.2.1 S u m m e r currents The water circulation arising from baroclinic pressure gradients computed from summer densityfieldsis discussed here. Figs 6.1 and 6.2 show the current velocityfieldobtained at 10-m depth for the June 24-July 7 1991 time period. Notice that model results will from now on be displayed for the zoom region, box shown on Fig. 6.1, only. The 40- and 100-m depth circulation patterns are respectively shown on Figs 6.3 and 6.4. The densityfields,from which the —* baroclinic forcing terms (Rres) were calculated, can be found on Figs 4.14 to 4.16 and the barotropic forcing term (R e*) on Fig. 4.17. The computed 10-m depth currents arising r from the other two summer densityfields,i.e. from the July 1-6 1985 and July 5-August 18 1991 time periods, are respectively shown on Figs 6.5 and 6.6. Two cyclonic gyres, covering the entire central-eastern part of Dixon Entrance, are among the interesting circulation features obtained for this season (see circles on Figs 6.2, 6.5, and 6.6 for the approximate location of the eddies). The combination of the two gyres would possibly correspond to the known Rose Spit Eddy (see Section 2.5). From the 1991 density data (Figs 6.2 and 6.6), two gyres of similar intensity were simulated while a strong eastern gyre and weak western one were obtained from the 1985 density data (Fig. 6.5). Note that the computed gyres are stronger in Fig. 6.6 than for the other two summer times. In order to estimate the rotation period of the Rose Spit Eddy, drifter simulations were performed at 10-m depth with the computed current velocityfieldcharacterizing the three summer periods. The 3D particle tracking algorithm written by Brian 0. Blanton 146 Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 134 W 134W 133W 133W 132W 132W 131 W 131 W 130W 130W Figure 6.1: Current velocityfield,at 10-m depth, driven by baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991. Units are cm s and the box delimits the zoom region shown on Fig. 6.2. _1 Chapter 6. Currents Arising from River Runoff and BarocUnic Effects 133 W 55 N 132 W 147 55 N - 54N -> -20 = 40 I 133 W 132W 131 W Figure 6.2: Current velocityfield,at 10-m depth, driven by baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991 for the zoom region (box shown on Fig. 6.1). Units are cm s and the two circles show the approximate position of the cyclonic eddies discussed in the text. -1 Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 148 133 W 55 N 54N 54 N 133 W 132W 131 W Figure 6.3: Current velocityfield,at 40-m depth, driven by baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991. The units are cm s . _1 Chapter 6. Currents Arising from River Runoff and Barochnic Effects 149 Figure 6.4: Current velocity field, at 100-m depth, driven by baroclinic forcing computed from the density field measured in the period June 24-July 7 1991. The units are cm s . -1 Chapter 6. 150 Currents Arising from River Runoff and Baroclinic Effects Figure 6.5: Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period July 1-6 1985. Units are cm s and the two circles show the approximate position of the cyclonic eddies discussed in the text. _1 Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 55 54 N- 151 •55 N -54 N N- 133W 132W 131 W Figure 6.6: Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period July 5-August 18 1991. Units are cm s and the two circles show the approximate position of the cyclonic eddies discussed in the text. _1 Chapter 6. Currents Arising from River Runoff and Baxochnic Effects 152 from University of North Carolina was used. A period of rotation of 6 to 32 days, 5 to 25 days, and 5 to 19 days was respectively obtained for the periods June 24-July 7 1991, July 1-6 1985, and July 5-August 18 1991. The variance in the rotation period, for a specific summer cruise, depends on the initial position of the numerical drifter (distance from the center of the gyres) and whether the drifter completes a loop around one of the two gyres or around a combination of the gyres. For example, using model results from June 24-July 7 1991, it takes on average 25 days to complete a full loop around the combined gyres and 6 to 10 days around only one of the two gyres. Using the 1985 model results, it takes 5 days for a numerical particle to travel around the strong eastern gyre and 25 days around a combination of the two gyres. Notice that observations suggest a rotation period similar to what obtained numerically here; that is 20-30 days (Crawford and Greisman, 1987) and 5-40 days (Crawford, 1995). In western Dixon Entrance, along the eastern side of Learmonth Bank, a strong northward current is calculated at all depths for the period June 24-July 7 1991. On the northern side of the bank, watersflowout of the entrance at the surface and into it at depth. This circulation pattern is however not found in the computed velocity field obtained with the other two summer cruises. A lack of density observations in 1985 (see Fig. 4.10) and the long time interval over which the density data were taken in the July 5-August 18 1991 time period could possibly contribute to this mismatch. Nevertheless, a common result to the three summer periods exist in western Dixon Entrance. The Learmonth Bank Gyre, observed in the tidal residual circulation, is completely absent from the three buoyancy flows. Chapter 6. Currents Arising from River Runoff and BarocHnic Effects 6.2.2 153 Fall currents Figs 6.7 and 6.8 display the 10-m depth buoyancy currents computed with the two sets of density observations available for the fall, i.e. taken respectively in the September 19October 10 1962 and October 14-21 1984 time periods. The two circulation patterns obtained for this season differ considerably from each others and from the summer time. A Rose Spit Eddy composed of three cyclonic gyres (see circles on Fig. 6.7) is simulated from the 1962 density field while no eddying motion is found in the 1984 computed currents. The main surface circulation pattern obtained for October 1984 I s a seaward (westward)flowover the entire entrance. Note that the breakdown of the Rose Spit Eddy, seen in model results computed for October 1984, w a s al- observed from current meter 30 observations (Bowman et al., 1992). In agreement with the summer results, no evidence of the Learmonth Bank Gyre is obtained from the fall cruise data. Drifter simulations were performed at 10-m depth using the model results obtained for the fall 1962. The Rose Spit Eddy was then found to have a rotation period lying in the 6 to 40 days range. This period of rotation is a bit longer than the ones obtained for the summer but still falls into the range estimated from observations. In 1967, Crean calculated, from density observations taken in September 19-October 10 1962, the distribution of dynamic height based on a 125-m reference level. There is considerable similarity between his result (Fig. 6.9) and the sea-surface elevation obtained with our model for the same time period (Fig. 6.10). 6.2.3 Spring currents The only density observations available for the spring time, that is in April 11 and 22 of 1984, have a poorly spatial coverage (see Fig. 4.12). As observations are totally absent from the western and eastern parts of Dixon Entrance, one can only have confidence in Chapter 6. Currents Arising from River Runoff and Barochnic Effects 55 N 154 55 N 54 N- 54 N 133 W 132 W 131 W Figure 6.7: Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period September 19-October 10 1962. Units are cm s and the three circles show the approximate position of the cyclonic eddies discussed in the text. _1 Chapter 6. Currents Arising from River Runoff and BarocHnic Effects 133 W 132 W 155 •55 N •54 N 133W 132 W 131 W Figure 6.8: Current velocityfield,at 10-m depth, driven by baroclinic forcing computed from the densityfieldmeasured in the period October 14-27 1984- The units are cm s . -1 Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 156 Figure 6.9: Distribution of dynamic height (m), between 0 and 125 m, for the period September 19-October 10 1962 (Fig. 28 in Crean, 1967). Contour interval is 0.01 m (1 cm). Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 133 W 132 W 157 131 W Figure 6.10: Sea-surface elevation generated by baroclinic forcing from the density field measured in the period September 19-October 10 1962. The thick solid line represents the zero surface elevation contour, thin solid lines positive values and dotted-dashed lines negative values. Contour interval is 1 cm. Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 158 55 N 54 N- -54 N 133 W 132 W 131 W Figure 6.11: Current velocity field, at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period April 11-22 1984- Units are cm s and the circle shows the approximate position of the cyclonic eddy discussed in the text. _1 model results of central D E . Fig. 6.11 shows the computed buoyancy currents for this time period. Only one cyclonic gyre (see circle on Fig. 6.11), located in the central part of the entrance, is present in the computed flow pattern. From numerical drifter tracking, a period of rotation of 13 to 48 days was found. This period is slightly longer than the one calculated for the summer and fall but a lack of observations might be the cause. Chapter 6. Currents Arising from River Runoff and BarocUnic Effects 6.2.4 159 Winter currents The 10-m depth buoyancy currents for December 2-18 1991, the only winter period for which density observations are available, are shown on Fig. 6.12. The Rose Spit Eddy is simulated here by a combination of two cyclonic gyres (see the circles on Fig. 6.12). Note however that those gyres are now confined to a smaller geographical region than during the other seasons. The rotation period of the eddy obtained numerically, at 10-m depth, varies between 3 and 27 days. As for June 24-July 7 1991, a net surface outflow, and deep inflow (not shown), are observed north of Learmonth Bank. There is also a strong outflow, present at all depths, characterizing the region north of Langara Island. Another feature is the strong westward currentflowingalong the Rose Spit sill. This current could have been set up by strong winter winds (see Fig 5.8). 6.3 Model Sensitivity To study its sensitivity to open boundary conditions, the model was rerun on the Fine Grid (Fig 4.3). It was found that moving the open boundaries closer to the Learmonth Bank region significantly alters, in the neighbourhood of the bank, the computed buoyancy currents. The model was also run on the Extended Grid (used for the results given in this chapter) with different boundary conditions (see Section 4.2); the depth-averaged flow was set to be tangential to all open boundaries except for boundary 5 (see Fig. 4.4) where it was set to be normal. Changing the direction of the depth-averaged flow at open boundaries has a negligible effect on the water circulation, only the solution near the open boundaries is affected. Therefore, the Extended Grid should be used over the Fine Grid and one should not be too concerned about setting the right direction of the depth-averagedflowat the open boundaries. Using a coarser grid did not modify the circulation pattern but revealed less detail. For Chapter 6. Currents Arising from River Runoff and Baroclinic Effects 54 N 160 54 N 133 W 132W 131 W Figure 6.12: Current velocityfield,at 10-m depth, driven by baroclinic forcing computed from the density field measured in the period December 2-18 1991. Units are cm s and the two circles show the approximate position of the cyclonic eddies discussed in the text. -1 Chapter 6. Currents Arising from River Runoff and BarocHnic Effects 161 example, for the period June 24-July 7 1991 and a depth of 10 m, only one cyclonic gyre was obtained in central-eastern Dixon Entrance when using the Coarse Grid (Fig. 4.2), compared to two with the finer grid. The vertical grid resolution affects the model results within the same order of magnitude as the friction parameters, that is about 5 cm s _1 away from open and land boundaries. Note that decreasing the friction parameter Cioo increases the depth-averaged currents but not necessarily the currents found at a specific depth. 6.4 Influence of Earth's Rotation and of Stratification on Eddy Viscosity The influence of the earth's rotation on the currents was obtained, as in the previous chapter for tidal residual and wind currents, by rerunning the model without rotation. An example of the numerical results is shown on Fig. 6.13. A large speed increase is one major difference between Figs 6.2 and 6.13. The maximum speed has actually doubled, away from open and land boundaries, from 30 cm s -1 on the rotating earth, to 60 cm s . _1 Note also that the circulation pattern is completely altered by the rotation of the earth. The influence of the stratification on the eddy viscosity was found to be, as for tides and wind, negligible on the buoyancy currents. 6.5 Summary The model simulations showed that : • The currents generated by river runoff forcing are small compared to tidal-rectified and buoyancy currents, except near the Skeena and Nass river mouths. • The buoyancy currents are Chapter 6. Currents Arising from River Runoff and Bcirochnic Effects 133 W •fc. K. N\ 132 W 162 131 W K. TV "N *\ ^ ^- , , t r r 11\v- * *• «r urrf"* it - *rfrf * 54 N 133 W 132W 131 W Figure 6.13: Current velocity field, at 10-m depth, driven by baroclinic forcing on a non-rotating earth. Units are cm s and the forcing term was calculated from the density field observed in the period June 24-July 7 1991. _1 Chapter 6. Currents Arising from River Runoff and Barochnic 163 Effects 1. sensitive to the location of open boundaries, but not to the depth-averaged conditions (flow direction) set at those boundaries, 2. not sensitive to the horizontal grid resolution, but more details are revealed on a finer grid, 3. as sensitive to the friction parameters as they are to the vertical grid resolution, 4. strongly influenced by the earth's rotation, 5. not influenced by the stratification effect on the eddy viscosity. • The Rose Spit Eddy (cyclonic gyre found in central-eastern Dixon Entrance) 1. arises from baroclinic effects, 2. is generally composed of two or three smaller cyclonic gyres, with seasonal variations in position and magnitude, 3. has a period of rotation, at 10-m depth, varying between 3 and 40 days, depending on the initial position of the numerical particle and the season (note that a short rotation period corresponds to a loop around one of the small gyres while a long one corresponds to a loop around a combination of the gyres), 4. is present in all seasons but can occasionally disappear, as in October 1984• In western Dixon Entrance 1. buoyancy currents show no presence of the Learmonth Bank Gyre for any of the seasons, 2. the computed flow pattern differs considerably from one baroclinic forcing (densityfield)to the other. Chapter 7 Evaluation of Model Accuracy The model results are expected to be a good approximation of reality. The purpose of this chapter is to evaluate the model ability to reproduce the real physical system. To have results as close as possible to observations, the four forcings (tides, wind, river runoff, and baroclinic pressure gradients) must be included in the model. The complete model solution, water circulation driven by all forcings, is first presented. Model results are next compared with two types of existing current observations. These are current data from moored buoys (current meters) and current data measured with surface drifters. From the former, we have temporal coverage and data for the we have spatial coverage and data for the surface deepflowwhile from the latter circulation. A summary concludes the chapter. 7.1 Mean Water Circulation Driven by all Forcings The model used in this thesis is weakly non-linear (see Section 3.5). That is to say the interaction term of the residual flow with itself is neglected. The model nonlinearity then lies exclusively in the tidal forcing terms and the complete model solution is simply a linear combination of the solutions obtained under the individual forcings. In other words, the full picture of water circulation is the sum of the water circulation patterns obtained in the previous two chapters. Note that forcing from tides, river runoff, and baroclinic pressure gradients were included in all model simulations presented in this chapter. However, for simplicity, the wind stress was only included for the 164 December Chapter 7. Evaluation of Model Accuracy 165 2-18 1991 time period; time when wind data were available and wind forcing was found to be non-negligible. It is nevertheless important to recall that wind effects are indirectly included in the density distributions used to compute the baroclinic forcing. Also, as the effect of stratification on the vertical eddy viscosity was found to be negligible with the four forcings, it was excluded from the calculations. 7.1.1 N u m e r i c a l results Numerical results for the June 24-July 7 1991 and October 14-27 1984 periods are shown on Figs 7.1 and 7.2, respectively. Note that results are displayed for the same ge- ographical region as the buoyancy currents found in the previous chapter (see box in Fig. 6.1). The current velocity fields shown in Figs 7.1 and 7.2 are mainly a linear addition of the tidal residual currents (Fig. 5.3) with the buoyancy currents (Figs 6.2 and 6.8 respectively). For all the time periods studied in this thesis, the known circulation features simulated are the Rose Spit Eddy and Learmonth Bank Gyre (see Fig. 7.1 for example). The only exception is the October 14-27 1984 period for which the Rose Spit Eddy is not observed and the Learmonth Bank Gyre is less well defined (see Fig. 7.2). From wind and current time series analysis, Bowman et al. (1992) also observed, in October 1992, a breakdown in the Rose Spit Eddy. They related this breakdown to the occurrence of the first winter storm in the region. It then seems that a strong fall wind event could considerably alter the gyral circulation normally present in Dixon Entrance. Drifter simulations performed numerically at 10-m depth showed that adding the tidal residual currents to the buoyancy currents generally increases the rotation period of the Rose Spit Eddy by a few days. Chapter 7. Evaluation of Model 133 W 166 Accuracy 132 W 131 W 55N 54N 132W 133 W 131 W Figure 7.1: Current velocity field, at 10-m depth., driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period June 24-July 7 1991. Units are cm s and the two circles show the approximate position of the Learmonth Bank Gyre (left) and the Rose Spit Eddy (right). _1 Chapter 7. Evaluation of Model 133 W 167 Accuracy 132 W 131 W - 132W 133 W 54 N 131 W Figure 7.2: Current velocity field, at 10-m depth, driven by tides, river runoff, and baroclinic forcing computed from the density field measured in the period October 14-27 1984- Units are cm s and the circle shows the approximate position of the Learmonth Bank Gyre. -1 Chapter 7.2 7. Evaluation of Model Accuracy 168 Comparison of Model Results with Current Measurements The evaluation of model accuracy consists in the comparison of the computed velocity field with the mean circulation observed from the physical system under study. This mean circulation can be brought to light by averagingfixed-pointcurrent meter data for periods longer than a few tidal cycles, or by Lagrangian tracer studies involving horizontal displacements of drifters over similar periods. When averaged over long periods, particle (Lagrangian) mean velocities are however likely to differ considerably from fixed point (Eulerian) means. Nevertheless, a persistent circulation of reasonably large amplitude would stand out in Eulerian as well as in Lagrangian observations. A list of the current measurements made through the years in the Dixon Entrance System can be found in Appendix B. To compare model results with observations, the model forcings and current observations should belong to a common time period. Therefore, any current observations taken during a time period covered by one of the seven CTD cruises considered in this thesis (see Figs 4.10 to 4.12) can be used to evaluate the model accuracy. 7.2.1 The 1962 current observations In the September 24-October 9 1962 time period, Ekman-type current meter measurements were taken, at five different depths with sampling interval of two hours, at three locations in central Dixon Entrance (Tabata, 1980) (Fig. 7.3). Note that the three anchor stations were occupied for only 2 to 5 tidal cycles each. Fig. 7.4 shows the average current components obtained from these measurements by Crean (1967). Model results, computed velocity profiles, for a similar time period (September 19-October 10 1962) and the same three locations are displayed on Fig. 7.5. Despite the short time interval over which the current measurements were taken, the computed current direction at both Chapter 7. Evaluation of Model Accuracy 169 Figure 7.3: Locations (A,B,C) where current measurements were taken in the period September-October 1962. Figure 7.4: Average current components at the three locations shown on Fig. 7.3 (Fig. 27 in Crean, 1967). Chapter 7. Evaluation of Model SPEED(cnVs) SPEED (cm/s) South/West -20 -16 -12 -8 -4 1 ^ — i 1 1 0 North/East 4 8 1 1 1 1 170 Accuracy South/West -20 -16 -12 -8 -4 1 1 I 1 1 1 0 SPEED(cm/s) South/West -8 - 4 . 0 4 North/East 4 8 12 /\ / \ l ' l ' j North/East 12 16 20 24 1 • 1 ' 1 1 1 1 /1 1 / 1 \• / / 1/ 8 1 "1 1 1 // '/ //// - n h - 1 i ( \ ) I Station A i i 1 • N i L L i i Station C 1 i i 1 i Figure 7.5: Computed velocity profiles at the three stations shown on Fig. 7.3. The model was driven by tides, river runoff, and baroclinic forcing arising from the density field measured in the period September 19-0ctober 10 1962. The dashed and solid lines respectively represent the west-east and south-north components. Chapter 7. Evaluation of Model 171 Accuracy stations A and C, respectively toward the south-west and north-east, are consistent with the observations. Note that station B was located near the expected center of the Spit Eddy, 7.2.2 Rose a less well defined current region. The 1984-85 current observations For the April-October 1984 a n < i October 1984-July 1985 time periods, 16 and 12 moorings were respectively deployed in Dixon Entrance and 1, in each case, in Hecate Strait. Three different densityfields,taken in 1984-85, were used in this thesis to drive the numerical model (see Figs 4.10 to 4.12). Note that they all correspond to either a deployment or recovery time of the moorings. Thus, the density and current observations do not completely overlap in time. To compare the current observations with the computed mean circulation, we first need to remove the high-frequency oscillations (including diurnal and semidiurnal tides) from the current time series. Note that our interest lies in low-frequency signals and thus any high-frequency components are considered as noise. Secondly, we need to average the current data over the time period covered by the model forcings, more specifically the density data. It should be pointed out that only the current meter times series which last more than 2 tidal cycles are considered here. A low-pass frequencyfilter,i.efilteringthe time series in the frequency domain, was chosen to attenuate the high-frequency waves. Its main advantage over a digitalfilterin the time domain is that all the high frequencies can easily be removed and the lower ones kept intact while no data are lost at the ends of records. Remember that density data exist only at the beginning and/or at the end of current meter records so that the number of data lost at the ends of the current time series should be minimized in order to retain some time overlap between density and current data. Notice that some care must be taken when applying this type offilter.First, the cutoff frequency should be chosen in a Chapter 7. Evaluation of Model 172 Accuracy band where no important energy lies. Secondly, a window taper should be used in order to minimize the Gibbs phenomenon present with sharp edges. Finally, at the ends of the time series, the signal amplitude should be as close as possible to zero. In order tofillall those requisites, the average and linear trend werefirstremoved from the time series, the discrete Fourier transforms calculated, a cosine taper applied such that all frequencies greater than 0.1 cycle per day and less than 0.05 cycle per day were respectively removed and retained, and the time series reconstructed from the frequencies left and by adding back the average and linear trend removed previously. Note that some verification was made to ascertain that ringing effects were negligible at the ends of time series. Figs 7.6 and 7.7 show the computed current vector field, at 50- and 150-m depth respectively, driven by tides, river runoff, and baroclinic forcing for April 11-22 1984- The observed average currents, at depths of 40-60 m and 140-160 m, are respectively displayed, with bold arrows at the center of a circle, on the same twofigures.Despite a temporal coverage for the current meters of only 50%, on average, of the cruise duration (26 to 69% and 33 to 67% at respectively 50- and 150-m depth), the model results and observations match reasonably well. One exception is the three (two at 150-m depth) current meters located on the western side of Dixon Entrance (arrows shown in a double circle on Figs 7.6 and 7.7). As this is a data-void region for density observations (see Fig. 4.12), no baroclinic forcing is available and thus the model is not expected to do well. The other exception is the current meter located right at the center of the entrance (arrow shown in a triple circle on Figs 7.6 and 7.7). The current direction is wrong at the two depths but not the speed. Fig. 7.8 shows the current comparison but for the October 14-27 1984 time period at 50-m depth. Here again the current observations cover on average only 50% (from 11 to 67%) of the cruise duration. Nevertheless, agreement of model results with observations is generally reasonable. Note that model results at 150-m depth, not shown here, also Chapter 7. Evaluation 133 W 133W of Model 173 Accuracy 132 W 131 W 132W 131W Figure 7.6: Current velocity field, at 50-m depth, driven by tides, river runoff, and baroclinic forcing computed from the density field measured in the period April 11-22 1984- The observed average currents are represented by bold arrows at the center of either a single, double, or triple circle and units are cm s . Note that current observations shown in double and triple circles are referred to in the text. _1 Chapter 7. Evaluation of Model 174 Accuracy Figure 7.7: Current velocity field, at 150-m depth, driven by tides, river runoff, and baroclinic forcing computed from the density field measured in the period April 11-22 1984- The observed average currents are represented by bold arrows at the center of either a single, double, or triple circle and units are cm s . Note that current observations shown in double and triple circles are referred to in the text. _1 Chapter 7. Evaluation of Model 175 Accuracy Figure 7.8: Current velocity field, at 50-m depth, driven by tides, river runoff, and baroclinic forcing computed from the densityfieldmeasured in the period October 14-27 1984. The observed average currents are represented by bold arrows at the center of either a single or double circle and units are cm s . Note that current observations shown in a double circle are referred to in the text. _1 Chapter 7. Evaluation of Model Accuracy 176 gives reasonable agreement with observations. I would like to emphasize the improvement of my model, for the Learmonth Bank region, over the Ballantyne et al. model (1996) (see Fig. 7.9). One of the main differences between Figs 7.8 and 7.9 is the circulation pattern characterizing the Learmonth Bank region. My model (arrows shown in a double circle on Fig. 7.8) gives a better match between numerical results and observations than does the Ballantyne et al. model (Fig. 7.9). Finally for the July 1-6 1985 time period, only few current meter observations were available. Comparison results are not shown here but a reasonable match between the model currents and observed ones was obtained. The only exception was near Learmonth Bank where there is a lack of density observations (see Fig. 4.10). 7.2.3 The 1991-92 current observations For the June-December 1991 and December 1991-July 1992 time periods, 5 and 4 moorings were respectively deployed in Chatham Sound, 3 and 1 in Hecate Strait, and 2 and 1 in Dixon Entrance. Heavy fishing restricted the deployment of further moorings in the entrance. The current time series were, as for the 1984-85 data, low-frequency filtered and averaged over time periods shared with density observations. Unfortunately, most of the current meters were installed in regions generally characterized by weak and spatially variable currents. Note also that the current meters in Chatham Sound, as in Hecate Strait, are expected to be strongly influenced by the winds (see Section 5.3). In the Dixon Entrance region, the few possible comparisons between numerical results and observations were generally satisfactory. As an alternative to current meters, 49 and 12 drifters were respectively deployed at 3- and 10-m depth during the summer 1991. From the 3-m drogued drifter positions, Ballantyne et al. (1996) calculated the average Eulerian currents (Fig. 7.10). Fig. 7.11 shows the computed current velocity field at the same depth (3 m) and time period Chapter 7. Evaluation of Model Accuracy 177 Figure 7.9: Velocities, at 50-m depth., computed by Ballantyne et al. (1996). Their model was forced by tides and baroclinic pressure gradients calculated from the October 14-27 1984 cruise. Observed currents are shown in bolder type with a small circle at their foot and the mooring name nearby. Single shafts in multi-shafted vectors represent 10 cm s . (Fig. 9(f) in Ballantyne et al., 1996) _1 Chapter 7. Evaluation of Model Accuracy 178 10 cm/s 25 cm/s Figure 7.10: Average currents as calculated from the summer 1991 drifter positions at 3-m depth. Estimated standard deviations are shown as ellipses around each vector. Lack of an ellipse indicates that only one observation was available (Fig. 10(a) in Ballantyne et al., 1996). Chapter 7. Evaluation of Model 133W 179 Accuracy 132 W 131 W 55 N -54 N 132W 133W 131 W Figure 7.11: Current velocity field, at 3-m depth, driven by tides, river runoff, and bare-clinic forcing computed from the densityfieldmeasured in the period July 5-August 18 1991. The units are cm s . - 1 (July 5-August 18 1991). There is a good match between Figs 7.10 and 7.11 except in Hecate Strait, where the flow varies strongly with the wind (Crawford et al., 1988). As a matter of interest, the model results obtained by Ballantyne et al. (1996), for similar time period and depth, are shown on Fig. 7.12. Note how much the two numerical results (Figs 7.11 and 7.12) differ in the western part of Dixon Entrance. From comparisons with observations (Fig. 7.10), we conclude here again that my model gives a better reproduction of the real physical system, especially in the Learmonth Bank region, than Chapter 7. Evaluation of Model 180 Accuracy Figure 7.12: Velocities, at 3-m depth, computed by Ballantyne et al. (1996). Their model was forced by tides and baroclinic pressure gradients calculated from the July 5-August 18 1991 cruise. Single shafts in multi-shafted vectors represent 10 cm s . (Fig. 10(b) in Ballantyne et al., 1996) _1 Chapter 7. Evaluation of Model Accuracy 181 does the Ballantyne et al. model. 7.3 Summary The model simulations, and comparison with current measurements, showed that : • The mean circulation of the Dixon Entrance System is primarily driven by baroclinic and tidal forcings. • The Rose Spit Eddy, arising from baroclinicity, is generally slowed down, by a few days, by tidal rectification. • The Learmonth Bank Gyre, arising from tidal rectification, is present in all seasons. • The model ability to reproduce the real physical system is high in Dixon Entrance (when the baroclinic forcing is known), undefined in Chatham Sound (not enough observations), and poor in northern Hecate Strait, where highly variable currents (due to wind effects) exist. • Fewer discrepancies between the physical system (or, more precisely, measurements of the physical system) and the model output are found with my model than with the Ballantyne et al. (1996) model, especially in western Dixon Entrance. Chapter 8 Conclusion Since as early as 1962, physical models have been used to study the physical oceanography characterizing the waters of the complex region that is northern British Columbia. A non-rotating barotropic, hydraulic, model was first developed. Then, barotropic finite difference models, either tidal and wind-driven, and a set offiniteelement models, a barotropic tidal model without and with a bottom boundary layer and a baroclinic diagnostic model, followed. This thesis constitutes an extension to the modelling effort started more than 35 years ago in order to better understand the dynamics behind the mean water circulation of the region. The model used here is an improved version of the baroclinic model described in Ballantyne et al. (1996). First, to adapt the model friction to the complex regional bathymetry, the parameterization of friction processes has been modified. Next, the stratification effect has been added to the vertical eddy viscosity formulation and the vertical advection of tidal momentum to the model. A refined and enlarged version of the horizontal triangular grid used by Ballantyne et al. (1996) has also been created. Finally, the wind stress and river runoff were added to the tidal and buoyancy forcing terms. The model ability to reproduce the real physical system has been proven to be higher with the model described in this thesis than its original version (the Ballantyne et al. model). Model results generally match the measurements of the physical system (Eulerian and Lagrangian current observations) well in Dixon Entrance, if density observations 182 Chapter 8. 183 Conclusion exist, but poorly in northern Hecate Strait. Note that highly variable currents are expected to be driven by the wind in Hecate Strait and a diagnostic model can then not well reproduce the currents of this region. Also, as too little information is available on the complex Chatham Sound region, it is not possible to determine how well the model is doing there. The model results were found to be • very sensitive to the horizontal grid resolution (bathymetric resolution) when driven by the tides, • sensitive to open boundaries when driven by the wind and baroclinic pressure gradients, • at least as sensitive to the bottom friction coefficient value (Cioo) as the vertical grid resolution, when driven by tides and baroclinic pressure gradients, • less sensitive to the coefficient for bottom and surface boundary layers (Ci) than Cioo• not sensitive at all to the inclusion of the stratification effect on the model friction. A grid made of the Extended Grid (Fig. 4.4) in the horizontal and 11-sigma levels in the vertical (Fig. 4.7(a)) is satisfactory for the purpose of this thesis. The best would be to use the Large Grid (Fig. 4.5) but the number of available density observations limits the extent of the grid we can actually use. Note however that to correctly model the wind-driven circulation, an inflow/outflow must be specified at the open boundaries of the Extended Grid. Thisflowwas here obtained byfirstrunning the model on the Grid. The physical insights provided by this study indicate that : Large Chapter 8. 184 Conclusion • The mean circulation of the Dixon Entrance System is primarily driven by baroclinic and tidal forcings. The average wind forcing is important in the winter time, but not in the summer; barotropic river runoff forcing is important only nearby the river mouths. • Important tidal rectification is generated, in order of importance, north of Langara Island, around Learmonth Bank, along Rose Spit Sill, and in the surrounding of Capes Chacon and Muzon. Topographic stretching/shrinking of tidal vorticity is the major vorticity generation mechanism along the sill, advection of tidal vorticity by the tides around the bank and Cape Muzon, while generation of vorticity by lateral variation in tidal vertical structure characterized the region north of Langara Island and around Cape Chacon. • Tidal rectification occurs along the Rose Spit Sill even on a non-rotating earth. • The vertical advection of tidal momentum by the tides decreases the tidal residual speeds to realistic values, mostly in Learmonth Bank region. Its inclusion in the model is essential. • The Rose Spit Eddy, a cyclonic gyre found in central-eastern Dixon Entrance, arises from baroclinic effects. Tidal rectification generally acts to slow it down by a few days. The eddy is commonly formed by the combination of two or three smaller cyclonic gyres and characterized by seasonal variations in position and magnitude. Its period of rotation, numerically obtained at 10-m depth, varies between 3 and 40 days, depending on the initial position of the particle (distance from the center of the gyres) and the season. Note that a short rotation period corresponds to a loop around one of the small gyres while a long one corresponds to a loop around a combination of the gyres. Chapter 8. 185 Conclusion • The Learmonth Bank Gyre, an anticyclonic gyre found in western Dixon Entrance, is generated by tidal rectification. The earth's rotation is essential for its formation. • Both the Rose Spit Eddy Rose Spit Eddy and Learmonth Bank Gyre persist in all seasons but the can occasionally disappear, such as during a fall weather transition in October 1984. • The residual currents are slowed down by the rotation of the earth and wind-driven current direction reversed in northern Dixon Entrance. The mean water circulation, or long-term drift, determines the movement of nutrients and pollutants, fish eggs and larvae, and disabled vessels. It is then very important to understand its dynamics, e.g. for the prediction of pollutant movement such as an oil spill. As the study region, northern BC waters, is almost entirely surrounded by land (islands and mainland), it constitutes an importantfisheriesand recreative area and the consequences of a disaster such as an oil spill could be considerable. 8.1 Future Work The model results could be improved through better model forcings and friction: • A spatially varying wind would be a better representation of the regional wind condition than the actual constant wind. • With the use of lighthouse wind data, the wind forcing could be included even for years previous to end of the 1980's, when ocean buoy data started. The land-based wind speeds, which do not well represent the ocean winds, would however need to be adjusted. A transfer function, estimated from land-based (lighthouse) and sea-based (ocean buoys) winds available for common time periods, could be used to adjust the data. Chapter 8. Conclusion 186 • The inclusion of river runoff input through baroclinic boundary conditions would better represent the reality than the actual barotropic conditions. • A better knowledge of the statistics of the density field (covariance function and mean) for the different seasons would improve the statistical interpolation process. Also, the inclusion of land boundaries in the interpolation process would be a major improvement. • Better knowledge of the sea floor sediments would allow the use of a spatially variable bottom friction coefficient, thus improving model friction and tided forcing. In brief, more information on the physical system, i.e. more observations, would definitively improve the model. Note also that the evaluation of the model accuracy could be improved with further current observations and model results with further grid refinement. I believe however that the future of this modelling effort lies in the application of a prognostic model, such as the finite element model developed by Lynch et al. (1996). The interaction of the residualflowwith itself and horizontal advection, assumed negligible in this study, could then be included. Also, as important internal tides are believed to exist in Dixon Entrance, it would be important to include them. A diagnostic model does not allow their inclusion. Furthermore, spring/neap tides are important in the study region and their associated increase/decrease in tidal turbulence could not be considered by our harmonic diagnostic model; the model assumes time-averaged tidal condition. Note that improving the model forcings would also be an asset for the prognostic model. Bibliography Aldridge, J.N., and A.M. Davies (1993). A high-resolution three-dimensional hydrodynamic tidal model of the Eastern Irish Sea. J. Phys. Oceanogr., 23:207-224. Atmospheric Environment Service (1981). Climatological Station Data Catalogue: British Columbia. Downsview, Ont., Canada. 94 pp. Atmospheric Environment Service (1982). Canadian Climate Normals: 1951-1980; Vol 5: Wind. Downsview, Ont., Canada. 283 pp. Bakun, A. (1973). Coastal upwelling indices, west coast of North America, 1946-71. NMFS-SSRF-671, NOAA. 103 pp. Ballantyne, V.A., M.G.G. Foreman, W.R. Crawford, and R. Jacques (1996). Threedimensional model simulations for the north coast of British Columbia. Contin. Shelf Res., 16:1655-1682. Barber, F.G. (1957). The effect of the prevailing winds on the inshore water masses of the Hecate Strait Region, B.C. J. Fish. Res. Bd. Can., 14:945-952. Barber, F.G., and A.W. Groll (1955). Current observations in Hecate Strait. Progr. Rep. Pac. Coast Stn 103, Fish. Res. Bd. Can. 23-25. Barber, F.G., and S. Tabata (1954). The Hecate oceanographic project. Progr. Rep. Pac. Coast Stn 101, Fish. Res. Bd. Can. 20-22. Bell, W.H. (1963a). Reproduction of estuarine structure and current observation techniques in the Hecate Model. Ms. Rep. Ser. 158, Fish. Res. Bd. Can. 50 pp. Bell, W.H. (1963b). Surface current studies in the Hecate Model. Ms. Rep. Ser. 159, Fish. Res. Bd. Can. 27 pp. BeU, W.H. and N. Boston (1962). The Hecate Model. Ms. Rep. Ser. 110, Fish. Res. Bd. Can. 35 pp. Bell, W.H. and N. Boston (1963). Tidal calibration of the Hecate Model. J. Fish. Res. Bd. Can., 20:1197-1212. Blake, R.A. (1991). The dependence of wind stress on wave height and wind speed. J. Geophys. Res., 96:20,531-20,545. 187 BIBLIOGRAPHY 188 Blumberg, A.F., and G.L. Mellor (1987). A description of a three-dimensional coastal ocean circulation model. In Heaps, N.S., editor, Three-Dimensional Coastal Ocean Models, volume 4 of Coastal and Estuarine Sciences, pages 1-16. AGU, Washington, D.C. Bornhold, B.D., and J.V. Barrie (1991). Surficial sediments on the western Canadian continental shelf. Contin. Shelf Res., 11:685-699. Bowden, K.F. (1978). Physical problems of the benthic boundary layer. Geophys. Surveys, 3:255-296. Bowden, K.F., and P. Hamilton (1975). Some experiments with a numerical model of circulation and mixing in a tided estuary. Estuarine Coastal Mar. Sci., 3:281-301. Bowden, K.F., L.A. Fairbairn, and P. Hughes (1959). The distribution of shearing stresses in a tidal current. Geophys. J. R. Astr. Soc, 2:288-305. Bowman, M.J., A.W. Visser, and W.R. Crawford (1992). The Rose Spit Eddy in Dixon Entrance: Evidence for its existence and underlying dynamics. Atmos.-Ocean, 30:7093. Bretherton, F.P., R.E. Davis, and C.B. Fandry (1976). A technique for objective analysis and design of oceanographic experiments applied to MODE-73. Deep-Sea Res., 23:559-582. Brower, W.A., H.F. Diaz, A.S. Prechtel, H.W. Searby, and J.L. Wise (1977). Climatic Atlas of the Outer Continental Shelf Waters and Coastal Regions of Alaska. NOAA, Boulder, CO. Vol 1: Gulf of Alaska, 437 pp. Celia, M.A., and W.G. Gray (1992). Numerical Methods for Differential Equations: Fundamental Concepts for Scientific and Engineering Applications. Prentice Hall, Englewood Cliffs, New Jersey. 436 pp. Cherniawsky, J.Y., and W.R. Crawford (1996). Comparison between weather buoy and comprehensive ocean-atmosphere data set wind data for the west coast of Canada. J. Geophys. Res., 101:18,377-18,389. Chevron Canada Resources Limited (1982). Initial Environment Evaluation for Renewed Petroleum Exploration in Hecate Strait and Queen Charlotte Sound. Crawford, W.R. (1996). Physical oceanography of the waters around the Queen Charlotte Islands. In Vermeer, K., R.W. Elner, and K.H. Morgan, editor, Ecology of Marine and Shoreline Birds of the Queen Charlotte Islands, British Columbia, in press. BIBLIOGRAPHY 189 Crawford, W.R., and P. Greisman (1987). Investigation of permanent eddies in Dixon Entrance, British Columbia. Contin. Shelf Res., 7:851-870. Crawford, W.R., and R.E. Thomson (1991). Physical oceanography of the western Canadian continental shelf. Contin. Shelf Res., 11:669-683. Crawford, W.R., A.V. Tyler, and R.E. Thomson (1990). A possible eddy retention mechanism for ichthyoplankton in Hecate Strait. Can. J. Fish. Aquat. Sci., 47:13561363. Crawford, W.R., M.J. Woodward, M.G.G. Foreman, and R.E. Thomson (1995). Oceanographic features of Hecate Strait and Queen Charlotte Sound in summer. Atmos.Ocean, 33:639-681. Crawford, W.R., W.S. Huggett, and M.J. Woodward (1988). Water transport through Hecate Strait, British Columbia. Atmos.-Ocean, 26:301-320. Crean, P.B. (1967). Physical oceanography of Dixon Entrance, British Columbia. Bull. Fish. Res. Bd. Can., 156:66 pp. Csanady, G.T. (1976). Mean circulation in shallow seas. J. Geophys. Res., 81:5389-5399. Cushman-Roisin, B. (1994). Introduction to Geophysical Fluid Dynamics. Prentice Hall, Englewood Cliffs, New Jersey. 320 pp. Daley, R. (1991). Atmospheric Data Analysis. Cambridge Univ. Press, Cambridge, USA. 457 pp. Davies, A.M. (1985). A three dimensional modal model of wind induced flow in a sea region. Prog. Oceanog., 15:71-128. Davies, A.M. (1986). A three-dimensional model of the north west European Continental Shelf, with application to the M 4 tide. J. Phys. Oceanogr., 16:797-813. Davies, A.M. (1990). On extracting tidal current profiles from vertically integrated twodimensional hydrodynamic models. J. Geophys. Res., 95:18,317-18,342. Davies, A.M. (1991). On using turbulence energy models to develop spectral viscosity models. Contin. Shelf Res., 11:1313-1353. Davies, A.M. (1993). A bottom boundary layer-resolving three-dimensional tidal model: A sensitivity study of eddy viscosity formulation. J. Phys. Oceanogr., 23:1437-1453. Davies, A.M., and J. Lawrence (1994). The response of the Irish Sea to boundary and wind forcing: Results from a three-dimensional hydrodynamic model. / . Geophys. Res., 99:22,665-22,687. BIBLIOGRAPHY 190 Davies, A.M., and J.E. Jones (1995). The influence of bottom and internal friction upon tidal currents: Taylor's problem in three dimensions. Contin. Shelf Res., 15:1251— 1285. Denman, K.L., and H.J. Freeland (1985). Correlation scales, objective mapping and a statistical test of geostrophy over the continental shelf. J. Mar. Res., 43:517-539. Dodimead, A.J. (1980). A general review of the oceanography of the Queen Charlotte Sound - Hecate Strait - Dixon Entrance Region. Technical Report 1574, Can. Ms. Rep. Fish, and Aquat. Sci. 248 pp. Donelan, M.A. (1982). The dependence of the aerodynamic drag coefficient on wave parameters. In Proceedings of the First International Conference on Meteorology and Air-Sea Interaction of the Coastal Zone, pages 381-387, Boston, Mass. American Meteorological Society. Flather, R.A. (1987). A tidal model of the northeast Pacific. Atmos.-Ocean, 25:22-45. Foreman, M.G.G., R.F. Henry, R.A. Walters, and V.A. Ballantyne (1993). A finite element model for tides and resonance along the north coast of British Columbia. J. Geophys. Res., 98:2509-2531. Fountain, A.G., and W.V. Tangborn (1985). The effect of glaciers on streamflow variations. Water Resour. Res., 21:579-586. Freeland, H.J. (1990). Sea surface temperatures along the coast of British Columbia: Regional evidence for a warming trend. Can. J. Fish. Aquat. Sci., 47:346-350. Gandin, L.S. (1965). Objective Analysis of Meteorological Fields. Israel Program for Scientific Translations, Jerusalem. 242 pp. Geernaert, G.L. (1988). Measurements of the angle between the wind vector and wind stress vector in the surface layer over the North Sea. J. Geophys. Res., 93:8215-8220. Geernaert, G.L., F. Hansen, M. Courtney, and T. Herbers (1993). Directional attributes of the ocean surface wind stress vector. J. Geophys. Res., 98:16,571-16,582. Geernaert, G.L., S.E. Larsen, and F. Hansen (1987). Measurements of the wind stress, heat flux, and turbulent intensity during storm conditions over the North Sea. / . Geophys. Res., 92:13,127-13,139. Gill, A.E. (1982). Atmosphere-Ocean Dynamics. Academic Press, San Diego. 662 pp. Glorioso, P.D., and A.M. Davies (1995). The influence of eddy viscosity formulation, bottom topography, and wind wave effects upon the circulation of a shallow bay. J. Phys. Oceanogr., 25:1243-1264. BIBLIOGRAPHY 191 Grant, W.D., and O.S. Madsen (1979). Combined wave and current interaction with a rough bottom. / . Geophys: Res., 84:1797-1808. Grenier, R.R., Jr., R.A. Luettich Jr., and J.J. Westerink (1995). A comparison of the nonlinear frictional characteristics of two-dimensional and three-dimensional models of a shallow tidal embayment. J. Geophys. Res., 100:13,719-13,735. Gross, T.F., and F.E. Werner (1994). Residual circulations due to bottom roughness variability under tidalflows.J. Phys. Oceanogr., 24:1494-1502. Haight, F.J. (1926). Nontidal currents in Southeastern Alaska. Geograph. Rev., 16:642646. Hamilton, K. (1984). Seasonal mean north Pacific sea level pressure charts, 1939-1982. Dept. Oceanogr., Univ. Brit. Col., Vancouver, B.C.,. Ms. Rep. 41: 177 pp. Haney, R.L. (1991). On the pressure gradient force over steep topography in sigma coordinate ocean models. J. Phys. Oceanogr., 21:610-619. Hannah, C G . (1992). The Influence of the Earth's Rotation on the Wind-Driven Flow in Hecate Strait, British Columbia. PhD thesis, University of British Columbia. 236 pp. Hannah, C.G., and D.G. Wright (1995). Depth dependent analytical and numerical solutions for wind-drivenflowin the coastal ocean. In Quantitative Skill Assessment for Coastal Ocean Models, volume 47 of Coastal and Estuarine Studies, pages 125- 152. AGU. Hannah, C.G., P.H. LeBlond, W.R. Crawford, and W.P. BudgeU (1991). Wind-driven depth-averaged circulation in Queen Charlotte Sound and Hecate Strait. Atmos.Ocean, 29:712-736. Harper, J.R. (1980). Coastal Processes on Norteastern Graham Island, Queen Charlotte Islands, B.C., volume Paper 80-1A, pages 13-18. Geol. Surv. Can. Hasselmann, S., et al. (WAMDI Group) (1988). The WAM Model: A third generation ocean wave prediction model. J. Phys. Oceanogr., 18:1775-1810. Henry, R.F., and R.A. Walters (1993). Geometrically based, automatic generator for irregular triangular networks. Commun. Numer. Methods Eng., 9:555-566. Hollister, H.J., and A.M. Sandnes (1972). Sea surface temperatures and salinities at shore stations on the British Columbia coast, 1914-1970. Pacific Marine Science report 72-13, Environment Canada, Victoria, B.C. BIBLIOGRAPHY 192 Hoos, L.M. (1975). The Skeena River estuary status of environmental knowledge to 1975. Special Estuary Series 3, Environment Canada. Hsu, S.A. (1981). Models for estimating offshore winds from onshore meteorological measurements. Bound.-Layer Meteor., 20:341-351. Hsu, S.A. (1986a). Correction of land-based wind data for offshore applications: A further evaluation. J. Phys. Oceanogr., 16:390-394. Hsu, S.A. (1986b). A mechanism for the increase of wind stress (drag) coefficient with wind speed over water surfaces: A parametric model. J. Phys. Oceanogr., 16:144150. Jamieson, G.S. (1985). Dungeness crab, Cancer Magisterfisheriesof British Columbia. Alaska Sea Grant Rep. 85-3, Alaska Sea Grant College Program, Fairbanks, Alaska. 37-60. Karl, T.R., and W.E. Riebsame (1989). The impact of decadalfluctuationsin mean precipitation and temperature on runoff: a sensitivity study over the United States. Clim. Change, 15:423-447. Kendrew, W.G., and D. Kerr (1955). The Climate of British Columbia and the Yukon Territory. Meteorol. Div., Dept. Transport, Toronto, 222 pp. Ketchen, K.S (1956). Factors influencing the survival of the lemon sole in Hecate Strait, British Columbia. / . Fish. Res. Bd. Can., 13:647-694. Kinney, P.J., J.C.H. Mungall, C.E. Abel, R.O. Reid, A.C. Vastano, and R.E. Whitaker (1976). Oil spill movement predictions. Technical Report 5, Kitimat Pipe Line, Ltd., Termpol submission re: marine terminal at Kitimat, B.C. Large, W.G., and S. Pond (1981). Open ocean momentumfluxmeasurements in moderate to strong winds. J. Phys. Oceanogr., 11:324-336. Le Traon, P.Y. (1990). A method for optimal analysis offieldswith spatially variable mean. / . Geophys. Res., 95:13,543-13,547. LeBlond, P.H., and L.A. Mysak (1978). Waves in the Ocean. Elsevier, Amsterdam. 602 pp. LeBlond, P.H., K. Dyck, K. Perry, and D. Cumming (1983). Runoff and precipitation time series for the coasts of British Columbia and Washington State. Ms. Rep. 39, Dept. of Oceanography, Univ. of British Columbia, Vancouver, B.C. BIBLIOGRAPHY 193 Liu, P.C., D.J. Schwab, and J.R. Bennett (1984). Comparison of a two-dimensional wave prediction model with synoptic measurements in Lake Michigan. J. Phys. Oceanogr., 14:1514-1518. Luternauer, J.L. (1976). Skeena Delta Sedimentation, British Columbia. Geol. Surv. Can. Paper 76-1A: 239-242. Ly, L.N. (1993). Effect of the angle between wind stress and wind velocity vectors on the aerodynamic drag coefficient at the air-sea interface. J. Phys. Oceanogr., 23:159-163. Lynch, D.R., and C E . Naimie (1993). The M tide and its residual on the outer banks of the Gulf of Maine. J. Phys. Oceanogr., 23:2222-2253. 2 Lynch, D.R., and W.G. Gray (1979). A wave equation model forfiniteelement tidal computations. Comput. Fluids, 7:207-228. Lynch, D.R., F.E. Werner, D.A. Greenberg, and J.W. Loder (1992). Diagnostic model for baroclinic, wind-driven and tidal circulation in shallow seas. Contin. Shelf Res., 12:37-64. Lynch, D.R., J.T.C. Ip, C E . Naimie, and F.E. Werner (1996). Comprehensive coastal circulation model with application to the Gulf of Maine. Contin. Shelf Res., 16:875906. Lynch, D.R.,and F.E. Werner (1987). Three-dimensional hydrodynamics onfiniteelements. Part LLinearized harmonic model. Int. J. Numer. Methods Fluids, 7:871-909. MacKay, B.S. (1954). Tidal current observations in Hecate Strait. / . Fish. Res. Bd. Can., 11:48-56. Mathews, W.H. (1958). Underwater gravel deposits of Dixon Entrance. Ms. Rept. 7, Univ. Brit. Col., Inst. Oceanog., Vancouver, B.C.,. 6 pp. McPhee, M.G., and J.D. Smith (1976). Measurements of the turbulent boundary layer under pack ice. J. Phys. Oceanogr., 6:696-711. Muccino, J . C , W.G. Gray, and M.G.G. Foreman (1994). Calculation of vertical velocity in a three-dimensional model using a least squares approach. Comput. Methods Water Resour., 10:1105-1112. Muccino, J . C , W.G. Gray, and M.G.G. Foreman (1996). Calculation of vertical velocity in three-dimensional, shallow water equation,finiteelement models. Int. J. Numer. Methods Fluids, submitted. BIBLIOGRAPHY 194 Munk, W.H. and E.R. Anderson (1948). Notes on a theory of the thermocline. J. Mar. Res., 7:276-295. Naimie, C.E., J.W. Loder, and D.R. Lynch (1994). Seasonal variation of the threedimensional residual circulation on Georges Bank. J. Geophys. Res., 99:15,96715,989. Ng, B. (1993). The prediction of nearshore wind-induced surface currents from wind velocities measured at nearby land stations. J. Phys. Oceanogr., 23:1609-1617. Pacanowski, R.C., and S.G.H. Philander (1981). Parameterization of vertical mixing in numerical models of tropical oceans. J. Phys. Oceanogr., 11:1443-1451. Pearson II, F. (1990). Map Projections: Theory and Applications. CRC Press, Florida. 372 pp. Pedlosky, J. (1987). Geophysical Fluid Dynamics. Springer-Verlag, New York, second edition. 710 pp. Petro-Canada (1983). Offshore Queen Charlotte Islands. Initial Environmental Evaluation. Calgary, Alberta. Phillips, N.A. (1957). A coordinate system having some special advantages for numerical forecasting. J. Meteor., 14:184-185. Pickard, G.L., and D.C. McLeod (1953). Seasonal variation of temperature and salinity of surface waters of the British Columbia coast. J. Fish. Res. Bd. Can., 10:125-145. Pond, S., and G.L. Pickard (1983). Introductory Dynamical Oceanography. Pergamon Press, New York, second edition. 329 pp. Powell, J. M. (1965). Annual and seasonal temperature and precipitation trends in British Columbia since 1890. CIR. 4296, CLI-34, Dept. of Transport, Meteorological Branch, Toronto,. 42 pp. Ricker, K.E., and J.W. McDonald (1992). Biophysical suitability of the North Coast and Queen Charlotte Islands regions of British Columbia for salmonid farming in net cages. Technical report, prepared for Province of British Columbia. Robinson, I.S. (1983). Tidally induced residual flows. In Johns, B., editor, Physical Oceanography of Coastal and Shelf Seas, chapter 7, pages 321-356. Elsevier, Amsterdam. 470 pp. Royer, T.C. (1979). On the effect of precipitation and runoff on coastal circulation in the Gulf of Alaska. J. Phys. Oceanogr., 9:555-563. BIBLIOGRAPHY 195 Royer, T.C. (1982). Coastal fresh water discharge in the northeast Pacific. J. Geophys. Res., 87:2017-2021. Ruddick, K.G., E. Deleersnijder, P.J. Luyten, and J. Ozer (1995). Haline stratification in the Rhine-Meuse freshwater plume: A three-dimensional model sensitivity analysis. Contin. Shelf Res., 15:1597-1630. Schwing, F.B. and J.O. Blanton (1984). The use of land and sea based wind data in a simple circulation model. J. Phys. Oceanogr., 14:193-197. Signell, R.P., R.C. Beardsley, H.C. Graber, and A. Capotondi (1990). Effect of wavecurrent interaction on wind-driven circulation in narrow, shallow embayments. J. Geophys. Res., 95:9671-9678. Simpson, J.H., and J.R. Hunter (1974). Fronts in the Irish Sea. Nature, 250:404-406. Slaymaker, O. (1990). Climate change and erosion processes in mountain regions of western Canada. Mountain Res. and Develop., 10:171-182. Smagorinsky, J. (1963). General circulation experiments with the primitive equations. I. the basic experiment. Mon. Wea. Rev., 91:99-164. Smith, J.D., and C E . Long (1976). The effect of turning in the bottom boundary layer on continental shelf sediment transport. Mem. Soc. Roy. Sci. Liege, t? serie, 10:369396. Smith, S.D (1980). Wind stress and heatfluxover the ocean in gale force winds. J. Phys. Oceanogr., 10:709-726. Smith, S.D (1988). Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature. J. Geophys. Res., 93:15,467-15,472. Snyder, R.L., M. Sidjabat, and J.H. Filloux (1979). A study of tides, setup and bottom friction in a shallow semi-enclosed basin. Part II: Tidal model and comparison with data. J. Phys. Oceanogr., 9:170-188. Soulsby, R.L. (1983). The bottom boundary layer of shelf seas. In Johns, B., editor, Physical Oceanography of Coastal and Shelf Seas, chapter 5, pages 189-266. Elsevier, Amsterdam. 470 pp. Stucchi, D.J., and U. Orr (1993). Circulation and water property study of Prince Rupert Harbour, summer 1992. Canadian technical report of hydrography and ocean sciences 154, Institute of Ocean Sciences. 42 pp. BIBLIOGRAPHY 196 Tabata, S. (1957). Heat exchange between sea and atmosphere along the northern British Columbia coast. Prog. Rep. 108, Fish. Res. Bd. Can. 18-20. Tabata, S. (1980). An inventory of physical oceanographic information for the waters of Queen Charlotte Sound, Hecate Strait, Dixon Entrance and their vicinity. PMSR 80-7, Inst. Ocean Sci., Dept. Fish. Oceans, Sidney, B.C. 24 pp. Tennekes, H. (1973). Similarity laws and scale relations in planetary boundary layers. In Workshop on Micrometeorology, pages 177-216, Boston, Mass. American Meteorological Society. Thompson, W.F., and R. Van Cleve (1936). Life history of the Pacific halibut. Report 9, International Pacific Halibut Commission. 184 pp. Thomson, R.E. (1981). Oceanography of the British Columbia Coast, volume 56 of Can. Spec. Publ. Fish Aquat. Sci. Canada Communication Group, Ottawa. 291 pp. Thomson, R.E. (1989). The Queen Charlotte Islands physical oceanography. In Scudder, G.G.E.. and N. Gessler, editor, The Outer Shores, pages 27-63. University of British Columbia. Thomson, R.E., and W.J. Emery (1986). The Haida Current. J. Geophys. Res., 91:845861. Trites, R.W. (1956). The oceanography of Chatham Sound, British Columbia. J. Fish. Res. Bd. Can., 13:385-434. Tyler, A.V. (1986). Hecate Strait Project: Results of the first two years of multispecies fisheries research. Technical Report 1470, Can. Tech. Rep. Fish. Aquat Sci. 50 pp. Visser, A.W., M.J. Bowman, and W.R. Crawford (1990). Dynamics of tidally forced basin-wide coastal eddies. In Cheng, R.T., editor, Residual Currents and Long-Term Transport, pages 64-78. Springer-Verlag, New York. Walters, R.A. (1992). A three-dimensional, finite element model for coastal and estuarine circulation. Contin. Shelf Res., 12:83-102. Water Survey of Canada (1988). Surface Water Data - Reference Index. Environment Canada, Inland Waters Directorate, Ottawa. 408 pp. Webster, I. and D.M. Farmer (1976). Analysis of salinity and temperature records taken at three lighthouse stations on the B.C. coast. Pac. Mar. Sci. Rep. 76-11, Inst. Ocean Sci., Sidney, B.C., Canada. 42 pp. BIBLIOGRAPHY 197 Westrheim, S.J., R.P. Foucher, W.R. Harling, and W. Shaw (1980). G.B. Reed groundfish cruise no 79-6, September 6-21, 1979. Technical Report 191, Can. Data Rep. of Fish, and Aquat. Sci. 64 pp. Xing, J., and A.M. Davies (1995). Application of three dimensional turbulence energy models to the determination of tidal mixing and currents in a shallow sea. Prog. Oceanog., 35:153-205. Xu, J.P., and L.D. Wright (1995). Tests of bed roughness models using field data from the Middle Atantic Bight. Contin. Shelf Res., 15:1409-1435. Zemba, J., and C A . Friehe (1987). The marine atmospheric boundary layer jet in the Coastal Ocean Dynamics Experiment. J. Geophys. Res., 92:1489-1496. Appendix A Glossary of Commonly Used Symbols Roman letters Ah, A Horizontal and vertical eddy viscosity coefficients A*, Al A in the bottom and the surface boundary layers for neutral stability A* Background eddy viscosity coefficient for neutral stability A® Total vertical eddy viscosity coefficient for neutral stability CJ,CD Bottom and surface (air-water interface) friction coefficients Cioo C) at 100 cm above the sea bottom CDN CD in neutral stratification C\ Coefficient used in the definition of A-, and A, / Coriolis parameter (= 2Qsin(f)) F(Ri) Reduction factor representing the influence of stratification on A° g Acceleration of gravity (= gz;g =9.8 ms ) G Barotropic pressure gradient vector (= —g^hVi) h(x,y), H Mean and total depth of the water (H — h + rj) h ,H Mean and total water depth of the model (h = h — z ) H = hm + vi) v v 9 -2 m m m Atmospheric forcing (— ^) v/^T i lb(z),l (z) Mixing length in bottom and surface boundary layers n Unit vector normal to the boundary AT Brunt-Vaisala. frequency B 198 0 m Appendix A. Glossary of Commonly Used Symbols p 199 Pressure —* Q Wave transport (from nonlinearity in continuity equation) R Baroclinic pressure gradient vector (= — ^- / ° VXpdz) Ri Richardson number t Time u(x,y,t) Horizontal water velocity vector with components (u,v) z —* U(x,y,z,t) Water velocity vector with components (u,v,w) Ub(x,y) Horizontal slip velocity near the seabed (100 cm above bottom) with components (ub,Vb). It is called bottom velocity. u\o{x,y) Horizontal wind velocity vector at 10 m height above sea surface u (x,y) Horizontal wind velocity vector with components (u ,v ) w w v%, ul w Friction velocity at the ocean bottom and surface (x,y,z) Unit vectors respectively to the east, north, and upward directions (x,y) Horizontal co-ordinates, positive eastward and northward respectively Xi,X Adjustable positive stratification parameters z Vertical co-ordinate, positive upwards with z = 0 at the mean sea surface 2 z Height of model bottom above the sea bed (=1 meter) Z Height above the sea surface ZQ Roughness length of the sea surface 0 Greek letters Aj,A, Thickness of bottom and surface boundary layers r)(x,y,t) Free surface elevation 7 Time-independent part of the bottom speed —* T Nonlinear part of the bottom stress /«o Von Karman's constant (= 0.4) Appendix A. Glossary of Commonly Used Symbols p, Coefficient of molecular viscosity v Kinematic molecular viscosity (= ^) uj Radian frequency of the n Jl Angular velocity vector of the Earth (= f2z; f2 = 7.292 x 10 <f> Latitude p Mass density of air th n a harmonic constituent p(x,y, z,t) Mass density of water po{z) Mean mass density of water p"(x,y,z,t) Perturbation mass density of water p» A reference mass density of water r , r , T" 6 s 1 200 -5 s ) Ocean bottom stress, ocean surface stress, and wind stress at ocean surface —* T Nonlinear advection term (Pj(x,y,z) Basis function associated with node j i?* Rotary components of u £o Length scale associated with the bottom boundary roughness c Amplitude of the free surface elevation at frequency (n) C Vorticity (= g - f ) (b Bottom vorticity (= fa - f*) n Special symbols V, V/i Gradient operator (J^, V , VI Laplace operator ( J ^ + J^- + J^-) and its horizontal component §- Total derivative (= §- + U • V) 2 t ', ( ) J^) and its horizontal component t t _1 Temporalfluctuationsand time average X Vertical average of X (= j-jU JH. Xdz) <> Domain integral (= J dxdy) hm & Appendix A. E ± Glossary of Commonly Used Symbols Rotary components (= (Ex ±iEy)/2) Appendix B Historical Review of Observations The purpose of this appendix is to review the measurements taken through the years in Dixon Entrance and surrounding waters. When modelling a real physical system, as it is the case in this thesis, any existing observation can be useful. Data related to the friction, forcings, and results of the model are considered here. Information on model friction can be acquired through seafloorsediments. Wind stress, river runoff, baroclinic pressure gradients, and tides are the important forcings. Corresponding observations are wind velocity data, river discharge rates, synoptic CTD observations from cruises, and sea level measurements through pressure gauges. Since the modelling of tides is not part of this thesis, the latter will be ignored here. Information on pressure gauges can however be found in Foreman et al. (1993). Finally, Eulerian and Lagrangian current measurements are used to evaluate model results. Other oceanographic observations, such as surface properties at lighthouses, temperature measurements from bathythermograph, wave measurements, and current measurements through ship-mounted Doppler current profiler have been made in the past. Only thefirstone of this category is discussed in this appendix. B.l Sea Floor Sediments Only a few data have been collected on sea bottom material of the region of interest. The earliest observations were taken by Mathews (1958) in Dixon Entrance. More recently, Westrheim et al. (1980) collected material in Hecate Strait and the Rose Spit area while 202 Appendix B. Historical Review of Observations 203 Harper (1980) conducted studies on surficial sediments along Mclntyre Bay (Fig. 1.1). Hoos (1975) and Luternauer (1976) also collected data but from the Skeena River Estuary. B.2 Winds Five types of wind observations are available for the north coast of British Columbia: 1. Synoptic Marine Observations The archive of synoptic marine observations constitutes a worldwide database to which volunteer ships-of-opportunity contribute (Petro-Canada, 1983). In this type of data, several biases could be present. For example, ships usually tend to avoid bad weather when possible and in that case data are biased toward good weather samples. Also, there are more observations in and near shipping lanes than in remote areas, leading to non-uniform sampling in space. Furthermore, the recorded winds vary with height of the instrument, as well as motion and speed of the vessel, and thus wind observations can easily differ from one ship to another. 2. Offshore Drilling Platforms An important negative aspect of wind observations taken from these platforms is the relatively small number of observations for the northern BC coast. In fact, the only wells that have been drilled were by Shell, none of them in Dixon Entrance, during April-November 1968 and March-May 1969 (Petro-Canada, 1983). Wind speed and direction were registered at the few sites but the time series are short. 3. Surface Pressure Maps Covering the 1945 to today period, six-hourly (monthly for 1945-1972 period) friction-modified geostrophic surface winds, calculated from surface pressure maps on a 3° grid are available (Bakun, 1973). Nevertheless, the only site close to the Appendix B. Historical Review 134 W of 204 Observations 133W 132 W •'i.iai .iiu-. <BJ. ji^ IB{|IJ - f' i 130 W 131 W Hiq"ii|ii!<l aP" dita T," ..HikLiCWlliV jijjil. , !i|; WSsi JIN Iijli!!|1jn' HS I , fly ' <f:jn5iai»irti|jb|li qi'lllliiii hi.Hl » JI I4:U!H|II-.„, IBI llftim. ««• Mii! 'KRiiiiii;; Tin 55 N- , •"•fmnnni*" "•yWh, .iiiiiitu l ' ' -S.ai.ir 'aiiilj. f. • Langara^. 54 N- T n p l 8 Rose.Spit .rlijlTfiiBnn. „,, Cpaen.ls-J Fit. BuoyM45 A liiJrBr'n S ill 1 -"nf .:l!l!»' B I> IE 4!HTH UI'I(INI <li>)lilf>«~i.ni|i1 jM'ililliH!! . .! t»i I i f i Masset ;||| ?j !j| ifjEj i!p :!j{jj|j|||!||i| iiiiiii ijniijj i i, i ,Pta, Lucy Is. Lucy Is H , '..'il&wjwla. *! P |, I i i"i'PiiMfirn !i | ( J " I| -54N i3ijamii r I. ' f":ii:!ifiHH. J J , !!t,t| ; "IsiSlallr ,^11111111111 H ill u, UI<1, lilll|8S»lf'' Buoy/183 - ' •llflljll ' ""L'lllllH. Shannon Bay " •tUf'tlohliolI fHLUn i!»..'-5' tttjanr U , r j. „,.. • .,i • „„„ ™±l Jn, c elS-M'iiltmmi uinuinliu • .•int l M , i | iiiililplilil 'ijjilinnlilj, j j|Mi|i]»i n f izr PartjgBmfi«s ,:a!lil i f v lSf '" i^i x !..!• "llfcil!!,!!!!,' , 'if ifeiijillllll lifeiW . If liH P* 3 g HI Oiimial •H J.i> 1 'ifili-i'n. .jiillliir - Buoy.205 Bar jffir'jiiii iiin: •lljli HlilBHililaB^13,1,11,!, 1 Bonilla.ls. taj,,™"!" , l'l'j 3 jl """ • i H, Sandspit ' l|i|iii| ::i .•'•!sai:Jliiii!5£d.lll.l,[|,a_ Ethelda Bay 53 H 134 W 133 W -53 N 132 W 131 W 130 W Figure B . l : Locations of land wind station (X), ocean buoys (filled triangle), and lighthouse stations ( • ) . study area, 54° N and 134° W , represents the oceanic winds and not the coastal winds. 4. Land Stations The Atmospheric Environment Service (AES) currently operates 8 climatological stations, and in the past operated an additional one, with anemometers for measuring wind speed and direction (Fig. B . l and Table B . l ) . Hourly wind records are provided. The U.S. National Weather Service also operates one station in the southern Alaska Panhandle. It is worthwhile to note that most of the land stations are here exposed to strong topographic effects caused by the orientation of Appendix B. Station Name Annette Island Bonilla Island Historical Review of Anemometer Period of Height (m) Record 1949-present Ethelda Bay Langara Island 13.7 Lawyer Island 47.0 Lucy 30.0 Island Prince Rupert Airport 13.7 or 10.0 Rose Spit 10.0 Sandspit (Automated) Triple Island 33.5 or 10.1 18.0 205 Observations Brief Site Description 1965-present Light beacon tower on rock NW of Bonilla Island. Good exposure in most directions, open water SE through S to N, 150 m hill 1km to E. 1957-present Hilly wooded island. Instrument on knoll 100 m SW of building, sheltered by hills and trees to E, S, and W. Terrain diverts N winds to NW and NE. 1967-present Lighthouse on NW tip of island, open sea from SW through N to NE. Dense woods NE through S to E. Instrument beside lighthouse. 1969-present Instrument on top of light beacon tower, good exposure in all directions. Everygreen trees 20 m high cover island. 1969-present Lighthouse on NE end of island. Good exposure from W through N to SSE. Wooded 30 m hill to SW shelters SW and W winds. Instrument west of terminal building, 1962-1966 surrounding area hilly and wooded. Open exposure to W, mountains to N, E and S. Small hill less than 1 km to E. 1971-present Flat low land, no vegetation. Excellent exposure in all directions. 1955-present Flat grassy area, exposed to water from NW through N to SE. Hills to W, SW. Thick trees to S and W. 1973-present On largest of three small rocks in Brown Passage, good exposure in all directions. Table B.l: Description of land wind station sites shown on Fig. B.l. At Rose Spit, an automatic wind sensing station is operated where wind direction and speed are radioed several times per hour to Prince Rupert. (From Atmospheric Environment service (1981; 1982) and Brower et al. (1977).) Appendix B. Historical Review of Observations 206 mountains and coast. 5. Ocean Buoys Three ocean data buoys are now maintained by AES in the surrounding of Dixon Entrance (Fig. B.l). These provide hourly oceanic winds at heights of 3.7 m and 4.9 m above water level. These data are however available only since late 1980'searly 1990's and the time series contain many gaps (see Cherniawsky and Crawford, 1996). In addition, two buoys with anemometers were deployed in Hecate Strait during the summer 1983 and winter 1983-84 by the Institute of Ocean Sciences and the Canadian Coast Guard (Crawford et al., 1988). Another anemometer was also deployed in April-October 1984 (Crawford and Greisman, 1987) . B.3 River Discharge River discharge observations have been made from Inland Waters Directorate of Water Survey of Canada since 1928 (Water Survey of Canada, 1988). Daily observations of freshwater flow rate are made by measuring the water level of the river, and using an empirical water level-discharge relationship specific to that river. The empirical relationship, characteristic to a specific river, is obtained by measuring the river width, water velocity, and thusflowratedischarge with corresponding water level. The relation between water level and discharge could change over the years and the calibration is thus redone once in a while. Discharge records of the Skeena and Nass rivers are maintained from gauge readings respectively at Usk (54° 37' 50" N; 128° 25' 40" W), 145 km upstream, and at Aiyansh (55° 12' 50" N; 129° 08' 20" W), 75 km upstream from the respective river mouth. Note that the locations of these gauges have been modified through the years. Both rivers Appendix B. Historical Review of Observations 207 have had a naturalflow(no upstream regulation) over the recording period, 1929-today for the Nass River and 1928-today for the Skeena River. Important gaps are present in both time series, that are 1949-1956 for the Nass River and 1932-1937 for the Skeena River. B.4 Oceanographic Cruises A summary of major cruises in the region of interest is presented in table B.2. To the best of my knowledge, the first oceanographic cruise dates back to September 1934; three stations were then visited by the University of Washington. Thefirstsubstantial oceanographic cruise, as well as the first Canadian cruise, was made in the summer of 1938 by the Pacific Oceanographic Group (POG). In 1948, the POG monitored the arrival and disappearance of the spring-summer run-off into and out of Chatham Sound. The purpose of the survey was to investigate the relationship between water properties and the migration pattern of salmon to their spawning grounds. Note that 1948 was a significant flood year, occurring once in 100 years for the Skeena and once in 25 years for the Nass, and hence the effect of full dilution was monitored throughout the Chatham Sound System (Ricker and McDonald, 1992). In 1954, a major oceanographic study, called the Hecate Project, was initiated by the POG. The goal was to define the water masses, their spatial and seasonal variability, and the tidal currents and circulation systems (Barber and Tabata, 1954). This was followed by the Coastal Surveys, in 1957, and the Coastal-Seaways Project, in 1958. The latter was undertaken to obtain data primarily to assess the changes in oceanographic conditions from the open ocean to the adjoining coastal seaways. In 1961, the Monitor Project, which was essentially a consolidation of the several previous projects into a single one, was undertaken. In 1967, the Oceanic-Coastal Program was initiated. This program was Appendix B. Historical Review Group* / Study* University of Washington* University of Washington* POG* POG* (Chatham Sound) POG* Hecate Project* Coastal Surveys* Coastal-Seaway Project* Monitor Project* Special Survey* of 208 Observations Time period September 1934 July 1937 May-June 1938 May-September 1948 May 1951 July-August 1951 May 1954 July 1954 August-September 1954 February 1955 April 1955 May-June 1955 November-December 1957 November-December 1958 April 1959 June 1959 November-December 1959 October 1960 February 1961 July-August 1961 September-October 1961 January 1962 March 1962 September-October 1962 Oceanic-Coastal Program* September 1967 NCODE* April 1984 October 1984 July 1985 PERD* June-July 1991 July-August 1991 December 1991 July 1992 Reference Dodimead (1980) Dodimead (1980) Dodimead (1980) Trites (1956) Dodimead (1980) Barber (1957) Crean (1967) Dodimead (1980) Dodimead (1980) Dodimead (1980) Dodimead (1980) Crean (1967) Dodimead (1980) Dodimead (1980) Thomson (1989) Ballantyne et al. (1996) Ballantyne et al. (1996) Not published Table B.2: Major oceanographic cruises in the Dixon Entrance System. Appendix B. Historical Review of Observations 209 designed primarily to provide information on the continental shelf and offshore waters off Vancouver Island northward to the QC Islands for two seasonal periods (early autumn and late winter). The data from this cruise were unique in that continuous profiles of temperature and salinity vs depth were obtained for the first time in the area. The POG, as it existed at the Nanaimo Biological Station, was slowly phased out, and in 1975 the Institute of Ocean Sciences (known as IOS) was established. Coastal Seaway work was resumed by IOS in 1977, and a new initiative called the Northern Coastal Oceanic Dynamics Experiment (NCODE) using modern computerized equipment was begun in 1983. The goal of NCODE was to study the tides and currents of northern BC waters to assist in the development of tidal numerical models and to determine environmental factors that may influence the groundfish and shellfish populations (Crawford et al., 1988). More recently, in 1991, the Panel of Energy Research and Development (PERD) has funded oceanographic studies in the area. The aim was to establish a basic understanding of surfaceflowsin the event of an oil spill; exploratory offshore drilling has indicated the possibility of oil reserves in the region. B.5 Current Measurements The first current observations made, to the best of my knowledge, were in 1931. At that time, the International Halibut Fisheries Commission used drift bottles to observe the surface waters (Thompson and Van Cleve, 1936). On four days during the summer 1948, Eulerian surface current observations, 10 to 40 hours long, as well as current observations at 4 intermediate depths were made in Chatham Sound (Trites, 1956). In Hecate Strait, thefirstcurrent observations to cover a complete tidal cycle were made in 1952; surface and bottom currents were then measured (MacKay, 1954). In the summer of 1954, hourly observations of currents, generally for periods of 50 hours, were made at various positions Appendix B. Historical Review Time period Summer 1948 July 1952 Summer 1954 Summer 1955 September-October 1962 May 1983-April 1984 April-October 1984 October 1984-July 1985 June-December 1991 December 1991-July 1992 of 210 Observations Number of Stations 3 in Chatham Sound 1 in Hecate Strait 1 in Dixon Entrance 6 in Hecate Strait 3 in Dixon Entrance 10 in Hecate Strait 16 in Dixon Entrance 1 in Hecate Strait 12 in Dixon Entrance 1 in Hecate Strait 5 in Chatham Sound 2 in Dixon Entrance 3 in Hecate Strait 4 in Chatham Sound 1 in Dixon Entrance 1 in Hecate Strait Reference Trites (1956) MacKay (1954) Barber and Groll (1955) Dodimead (1980) Crean (1967) Crawford et al. (1988) Crawford and Greisman (1987) Crawford et al. (1988) Bowman et al. (1992) Ballantyne et al. (1996) Ballantyne et al. (1996) Table B.3: Eulerian current measurements in the Dixon Entrance System. at the surface, bottom, and many intermediate depths (Barber, 1957). The first program of current meter observations devoted to obtaining a clear picture of the regional circulation was carried out in 1982-1985 (Crawford et al., 1988). Before that time the bulk of the current meter studies concentrated on measuring tidal currents for navigation purposes. A summary of Eulerian current measurements in the region of interest is presented in Table B.3. In the last decade, a number of Lagrangian measurements have been made in the study area. These include a satellite-tracked drifter, drogued at 10 m, followed during October-December 1983 (Crawford and Greisman, 1987). In order to study near-surface currents in Dixon Entrance, 11 Loran-C retransmitting drifters, drogued at 10 m, were launched in June 1984 (Crawford and Greisman, 1987) and 17 in September-October 1986 (Bowman et al., 1992). In the winter 1991, another satellite-tracked drifter passed in the Appendix B. Historical Review of Station Name Bonilla Island Green Island Langara Island Masset Port Clements Prince Rupert Sandspit Shannon Bay Triple Island Observations 211 Period of observation April 1960-present February 1935-August 1936 November 1936-August 1937 March 1940-present January 1940-October 1942 October 1941-August 1942 January 1940-May 1942 August 1953-December 1956 January 1940-August 1941 November 1939-December 1970 Table B.4: Lighthouse stations (see Fig. B.l for the geographical locations) where SSS and SST are, or have been in the past, measured with their period of observations. study area (Hannah, 1992). Under PERD, 49 and 12 drifters, respectively drogued at 3-m and 10-m depth, were deployed during the summer 1991 (Hannah, 1992; Ballantyne et al., 1996 ). In the winter 1992, 3 additional drifters were deployed, under PERD, in the surrounding of Dixon Entrance. Finally in the last few years, as part of the World Ocean Circulation Experiment (WOCE), some drifters passed near the region of interest. B.6 Lighthouse Observations Since 1934, daily observations of sea surface sahnity (SSS) and temperature (SST) have been made at various locations along the northern British Columbia coast (Fig. B.l and Table B.4). The observations, taken at a depth of 0.9 m, have routinely been made within 1 hour before the time of high tide. Note that over the recording years, measuring procedures have changed (Hollister and Sandnes, 1972). Appendix C Governing Equations The purpose of this appendix is to present the continuity equation and 3D shallow-water equations of motion for an incompressible and Boussinesqfluid,and briefly describe how they are derived. Pedlosky (1987), Gill (1982), Pond and Pickard (1983), CushmanRoisin (1994), and LeBlond and Mysak (1978) may be consulted for a more detailed derivation. C.l Continuity Equation for an Incompressible Fluid In the absence of sources and sinks of mass within afluid,the condition of mass conservation is expressed by the continuity equation : ^ + v . ( i o = o, P (ci) where p is the mass density of water, t the time, and U the water velocity vector with u,v,w components to the east, north, and upward directions respectively. (C.l) may also be written ?£ + (V.U) p = 0, (C.2) where ^ = §i-\-U • V, i.e. local time rate of change plus advection. For an incompressible fluid the conservation of mass becomes conservation of volume, the sound waves are then eliminated and (C.2) becomes the reduced continuity equation: V •£/ = (). 212 (C.3) Appendix C. Governing Equations C.2 213 Navier-Stokes Equations The equations of motion are the expression of Newton's second law of motion stated per unit volume for a fluid continuum : ~8U +v-(uu) - V p + pV$ + F, dt (C.4) where p is the pressure and $ the potential by which conservative body forces such as —* gravity can be represented. F is any nonconservative force such as the frictional force in a fluid. For a Newtonianfluidlike water, viscous stresses are proportional to velocity gradients (linear stress-strain relation) and F can be expressed as (Pedlosky, 1987) F = /zV <7-r^V(V-<7), 2 o where V = 2 (C.5) + J^- + J-^, i.e. the Laplace operator, and p is the coefficient of molecular viscosity. Equation (C.4) is stated to be in conservative form. A non-conservative form can be obtained by applying the continuity equation (C.3) to (C.4): Now introduce (C.3) into (C.5) and substitute the latter in (C.6). The Navier-Stokes equations of motion for an incompressiblefluid,i.e. —* ^ = - - V ] - g + uV U, 2 P (C.7) are thus obtained. g[= gi = V$;# =9.8 ms ) is the gravitational acceleration and v -2 (= ^) the kinematic molecular viscosity. C.3 Reynolds Equations The Navier-Stokes equations describe the motions of viscousfluidsfor all scales; from very large scales, such as ocean currents, to very small scales, such as capillary waves. Appendix C. Governing 214 Equations The turbulent motions transfer energy from the larger to the smaller scales where the kinetic energy is ultimately dissipated by molecular viscosity. The transfer of energy comes from the non-linear terms of the equations where all scales interact. Obviously the description in detail of all scales of motion is not realistic. Following Reynolds' procedure, equations for motions averaged over a suitable time period are considered. The velocity, pressure, and density fields are split into a mean (( )t) and a fluctuating or turbulent ( ') components, such as U — (U) + U', p = {p)t + p ' , and t P ~ (p)t ~r p'• B y definition, (U')t — 0. Substituting those components into (C.7) and taking the time average gives —r ^DT = ~ w v { p ) t t - ° + u V 2 { t I ) t - { & • v ( f , ) t - ( c - 8 ) That is the Reynolds equations. The last term on the right hand side of (C.8) represents the effect of the turbulent small scales on the mean flow. Assuming an incompressible fluid, (U'ffi-U')) t = 0 can be added to the turbulence terms of (C.8) without changing the value, only the mathematical form. For the x-component, the turbulence terms become d(u'u% d(u'v') d(u'w') dx dy dz t t These terms are the gradients of Reynolds stresses. Turbulent motions act on the large-scale flow in a manner that mimics the way in which molecular motions affect the macroscopic flow. Reynolds stresses can then be related to the mean flow gradient most simply by some turbulent or eddy viscosity coefficients as, for example, -<.V>, = -<*V), Ox = A J $ & ; Oy -<»V) = 4 , ^ , ( (CIO) Oz where the coefficients Ah and A are respectively called the horizontal and vertical turv bulent or eddy viscosity coefficients. 215 Appendix C. Governing Equations From now, only the mean equations of motion will be considered and for simplicity the brackets with subscript t will be omitted. Introduce (C.10) into (C.9), and thereafter into (C.8). The resulting Reynolds equations are — where VX = =--Vp- 9 +v• ( w ) + h T m [A.-J , (C.ll) + J^y, i.e. the horizontal gradient operator. Note that values for eddy viscosity are up to 10 times those for kinematic molecular viscosity (Pond and Pickard, 11 1983) and therefore the latter is neglected. C.4 Reynolds Equations in the Earth's Rotating Frame In a rotating frame such as the Earth, (C.ll) becomes ^ + 2(ixU = -^Vp-g + V -{A V U) h h h +^ [^J^j ' ( - ) C 12 where €l is the earth's angular velocity vector (= fi£;fi = 7.292 x 10 s ). -5 -1 If the horizontal area being considered on the earth is not too large, a tangent plane to the earth can be used; this is called the f-plane approximation. On a f-plane, (C.12) becomes Du si Dv , } v 1 dp = , „ ' + V k I dp . A „ , { A k V h U ) + „ .. „ . d I A du\ & M d ( .„ „„. • dv\ m -i * - + & j• -bl - - i ~ * • " + ai [ T,) • Dw +fu = + ldp p „ 9+ v . iA Vkw) A d (, < C 1 3 ) (C14) dw\ A (C16) P The Coriolis parameter, f = 2Vlsin</>, is the component of the planetary vorticity normal to the earth's surface at latitude <f> and is here taken to be constant and equal to its value at the center of the study area. Appendix C.5 C. Governing 216 Equations Boussinesq Fluid Introduce a perturbation density p", such that p — po{z) + p"(x,y,z,t). The density variations in the world ocean are relatively small, usually in the order of 0.1% the average density p (LeBlond and Mysak, 1978). Now assume there is some constant reference 0 value, p„, from which the fluid density does not depart much, p, could be thought as a depth-averaged value of po(z). In other words, the variations caused by the existing stratification and/or thefluidmotion are small compared to the reference value. In a first approximation, variations influiddensity can be neglected where they affect the inertia terms, but retained in the buoyancy terms where they provide forcing for density-driven flows. That is called the Boussinesq approximation. Note that it is a good approximation forflowsencountered in fresh and salt water but not forflowswith heavy sediment loads, where variation in density is of the same order as density. For a Boussinesqfluid,(C.13) and (C.14) respectively become Du . Di ~ f Du Dt v 1 = ~7.£ f u _ = -7.i . . _ + * • v 1 dp . + dp _ + v . d ( . du\ ei K + ., _ " - « ^ » , ; ) • ._. ., . < C 1 6 > d ( . dv\ + T> {^ei)' ( C 1 7 ) and (C.15) remains unchanged. C.6 Shallow-Water Equations Assume that motions of interest are those with a horizontal scale large compared with the vertical scale, i.e. shallow-water approximation. In this context, the vertical acceleration of the fluid is so small compared to gravitational acceleration that the Archimedian principle for a static fluid can be applied; that is the pressure difference between any two points on the same vertical line depends only on the weight of the fluid between Appendix C. Governing 217 Equations those points. The weight of thefluidis thus assumed to identically balance the vertical pressure gradient, that is the hydrostatic approximation, and (C.15) can be written as §f = - » . (C.18) The integration of (C.18) with respect to z gives p(z) = p(rj) + gp.n + g pdz, (C.19) where 77 is the free surface elevation and p(rj) the atmospheric pressure. ^(77) is here assumed to be constant. Now differentiating (C.19) with respect to x and y and introducing the result into the horizontal momentum equations (C.16-C.17), we obtain the shallow-water equations of motion: Du ~Dt Dv ~Dt Appendix D Derivation of the Model Equations The model equations are derived by applying diverse approximations and mathematical transformations to the governing equations presented in appendix C. The mathematical methods employed here are those introduced by Lynch and Werner (1987) and thereafter successfully used in Walters (1992), Lynch et al. (1992), Lynch and Naimie (1993), and Naimie et al. (1994). Apart from the friction terms, the model equations are those presented in Walters (1992) and are, for completeness, re-derived in this appendix. D.l Representation of Mixing Processes Kinetic energy is transferred, via turbulent motions, from larger to smaller scales. Large eddies break down into successively smaller eddies until the eddies are small enough for molecular viscosity to take effect. The energy is then dissipated into heat. In the ocean, energy is dissipated into heat at typical eddy sizes of the order of a millimeter (Gill, 1982). Therefore numerical models face a problem: it is not possible to solve all scales down to the millimeter. A common technique to overcome this problem is to make the viscosity artificially large (i.e. introduce an eddy viscosity) so that sufficient energy dissipation can occur on scales resolvable by the numerical scheme. The purpose of introducing eddy viscosities is then to simulate the physically real net cascade of energy from the larger than grid-size scale to the smaller scales which are truncated by the discrete differencing of the model. With the assumption that the grid scale lies within an inertial sub-range, i.e. there is 218 Appendix D. Derivation of the Model Equations 219 a net transfer of energy to higher wave numbers in the neighbourhood of the grid scale, horizontal diffusivity (Ah) can be expressed in the form (Smagorinsky, 1963): A = h (k S) : h \( du dv\ (dv fc-fljj {ai + du\ + „. T )' ( y D ' 1 ) where kh is an horizontal diffusion constant (~ 0.1 — 1.0) and 5 the grid size. Ah then depends on local horizontal shear and grid scale. As the grid scale decreases, Ah diminishes, and the need for horizontal viscous terms reduces. In the present study, the model grid is assumed sufficiently fine in the horizontal that horizontal diffusion terms are assumed negligible. Ah is here set to zero. Thus, unresolved mixing processes are from now on represented as vertical diffusion only. Using this approximation, the shallow water equations of motion (C.20-C.21) can be written, in conservative form, as du „ d(uw) ? _ „ d ( „ du\ g f° _ , — + V . {uu ) + -± - + fxu + gVhTi\ »ir )=-/ Vhpdz, dt dz dz \ dz J p Jz A J fc r ir „. (D.2) t where u = (u,v), i.e. the horizontal water velocity vector. Neglecting for now the stratification effect, the vertical diffusivity is expressed, using (3.31), (3.24) and (3.28), as : A = A + A' = «o(C o)^6(z)|<4| + « {j DNJ b v c v v 10 0 2 h(z)\u \. w (D.3) The continuity equation (C.3) remains the same : V -u+^ h D.2 = 0. (D.4) Bottom and Surface Boundary Conditions In order to solve the governing equations (D.2) and (D.4), it is necessary to apply appropriate boundary conditions. At the model bottom (100 cm above the seabed), the 220 Appendix D. Derivation of the Model Equations kinematic and stress boundary conditions are respectively given by w ( dh dh \ m -[ -d - = U m -dy-) + V X ( D 5 ) anc - T OZ -, (D.6) p, where h is the mean water depth of the model (h = h — 1 m) and f* the bottom stress. m m Using (3.11), equation (D.6) can be rewritten as A QV (D.7) C \u\u. = Z wo Note that from now on, the horizontal velocity at the bottom of the model (z = —h ) m will be called bottom velocity and be written as UbAt the free ocean surface, the same conditions are respectively given by and p» OZ where r*" is the wind stress. Using (3.14), we can rewrite (D.9) as A^ v oz = ^C \u \u p, DN w w = H $, m (D.10) where H (= h + 77) is the total water depth of the model. For convenience, a new m m notation for atmospheric forcing (H ^) was here introduced. m D.3 Vertically Integrated Equations The primitive continuity equation and horizontal momentum equations are vertically integrated from the model bottom z = —h (x,y) to the sea surface z = rj(x,y,t). m Appendix D. Derivation of the Model Equations 221 The vertical integration of the continuity equation (D.4), using (D.5) and (D.8), gives 5!Uv ot f c . r udz = ^ + V -(h J-h at h m + r,)Z=0. (D.ll) m Note that an overbar indicates the vertical average of a quantity, defined as the integral of a quantity from bottom of the model domain to the free surface divided by the total depth H , i.e. m — X = I f I — Xdz. rim •*—hm The vertical integration of the horizontal momentum equation (D.2), using (D.7) and (D.10), gives du ———— d(uw) -* „ Cioo\ub\ub ? "5 — + V • (uu) + + / xu + gVhtj + ' *" = * + R, 100 h Ot Oz H ,^ (D.12) m where a shorthand notation has been introduced for the baroclinic pressure gradient: R=—?-[ - r V dz. (D.13) hP P* Jz p D.4 Frequency Domain The free surface elevation and the three components of velocity can be expressed in terms of harmonic expansions in time : *(*.v.*) = 5 £ Z U(x,y,z,t)=l (D.14) Un(x,y,z)e- ^ , (D.15) iUnt n=-N £ Z <;n(x,y)e- , i t n=-N where N is the number of harmonics selected, w the radian frequency of the n harth n monic constituent with complex amplitudes (? , U ), and i — n n —* It is assumed that —* u;_ = — uj and correspondingly (c_ , C/_ )=(c , f/ )*, with the asterisk indicating comn n n n n n plex conjugation. It is also assumed that n = 0 represents the time-independent constituent, with LOQ = 0 and 5(<j, Uo) — 0. ^s(x) represents the imaginary part of x. 0 Appendix D. Derivation of the Model Equations 222 Equivalent expressions to (D.14-D.15) are r]{x,y,t) = s ,(x,y) + » (jT &.(*,y)e-' ^ , w re U(x,y, z,t) = U (x,y,z) + M (f) y, z ) - ^ ret , e where 3t(x) is the real part of x. The residual component is here defined as ( ^ . , £ U ) = ^(ft,tfo). (D.16) Note that the baroclinic forcing (R) and the wind stress (H $) can also be expressed in m terms of harmonic expansions in time : R = \ Z E He-'""', n=-N In this thesis, only the residual component (n = 0) of both R and i f ^ are considered. m The previous expressions then become R= -R\ l H* m D.4.1 (D.17) = R\ „ e = ^p*„ = H % . m (D.18) et Absolute value term Substituting (D.14-D.18) into the governing equations is straightforward for the linear and bilinear terms. For the absolute value term of the bottom velocity, the mean and time-variable parts can be separated as in Snyder et al. (1979): 1^61 = {itb-Ub)* (D.19) . *p=-N *p=-N q=-N Appendix D. Derivation of the Model Equations 223 where from now on the subscript 6 is written for convenience as a superscript. Using the Taylor series expansion, (D.19) can be truncated beyond quadratic terms and re-written as \u | = 7 (D.20) ° ' p=-N q=-N where 7 is here defined as the time-independent part of the bottom speed defined as the root-mean square component of bottom velocity : 1 N 7= D.4.2 > * , . • Ure, + 2 IX-" -, (D.21) 6 e G o v e r n i n g equations The introduction of harmonic expansion into the governing equations (D.2), (D.4), (D.ll), and (D.12) gives the following set of equations : 1 N dw dz (D.22) n n=-N i 1 1 - Y, [-w»«n + V -(fc iZ )]e- -»« = - V N i fc ra fl f c (D.23) .^ n=-N N -iw u n n d_ ( idiln + f xu + gV q n h dz\ n v -iu t _ n dz Ro + T + dT n n (D.24) 1 2 Z N S • ~=i , £ -=> , T7 , ^1007 -* H„ * 0 + #o + f T»=-JV n (D.25) where thefirst-ordervertical viscosity is Al = Ko(Cioo)'4(*)7 + «o I ^°\ Hm ) l,(z), (D.26) nonhnear forcing terms 1 N N = -9 £ Z £ ,7 ep->( "j>+'"«) Su t P p=-Nq=-N q t (D.27) 224 Appendix D. Derivation of the Model Equations N T» = - i N £ d(*w) £ ,-t(a/p+w,)t (D.28) ^ p=-Nq=-N L and high-order friction terms (D.29) p=-N q=-N N N E z_/ E N '100 g-, E °' p=-N q=-N a=-N >=-N ^ • g/ ^ x( P e " i K + a , , + w (D.30) ' ) t 9^-P The stress boundary conditions (D.7) and (D.10) are now given by I E K ^ e - ^ = \ [HA - t) {» = r,), (D.31) n=-N 1 N 2 S Z A/ — ^(f -f ) C7ioo7«n n (D.32) (z = -h ), n m az n=-JV L the kinematic boundary conditions (D.5) and (D.8) by (D.33) p=-Nq=-N n=-iV (z = Z n=-iV Z -h ), m (D.34) n=-JV and the boundary condition (3.10) by j If - E 2 D.5 • e '""'= 0 - — land boundary). ({x,y) (D.35) n=-7V Steady Component In this present work, the phenomenon of interest is the residual circulation and thus a steady state solution is sought. For the steady component, (D.22-D.35) are given by Vfc • U . + re V fc • {H u ) m rea dw re (D.36) 0, d2 = V •Q h r e s , (D.37) Appendix D. Derivation /x u + rea f xu of the Model Equations - j- ( A I ^ J V G , H R E A + gV <; rea h rea z 225 = R\ , + E f , + r e + % ^ t £ . = * . + Rres e Pe rim z (D.38) + f , - ^p, rim (D.39) r e where the first-order vertical viscosity is Al = Ko(C o)"/6(*)7 + "0 (H \$re,\) h(z), k 10 m (D.40) the nonlinear forcing terms Qr.s = ~ \ £ ^ - p , (D.41) p/0 I ^ - i E j v W v , ) ^ ) , (D.42) _ «o(Cioo)^ , , , ^ ^ ^ ^fltt-p T , = — k{z) 2s i • ™'>~oT> (UAd) P p#0 the high-order friction terms u Pe m u P p= —N > ^ p^O = ^4 E « ' < « ) ^ (D.44) t pfr^ p^O and the boundary conditions A l ^ +T Oz rea =H* m rea (z = r ), (D.45) (z = -h ), (DAG) (.Z = 7?), (D.47) (z = ~h ), (D.48) d boundary). (D.49) } r\ —* 4 ^ + f , = C7 07«re. + f , Pe 10 m Pe OZ 1 N Wre. = (Ure. ' ^hSre.) + J £ (u • V ^ _ ) p p p^O We* = ~ (Ure. ' Vhh ) P u -n rea m =0 m (( ,y) = ' x a r | Appendix D. Derivation of the Model Equations 226 Note that we are considering only a weakly nonlinearity here. That means the residual component is considered sufficiently weak that its influence on the nonlinear advection terms is neglected. Nevertheless, all continuity and frictional nonlinearities involving the residual component are included. D.6 Rotary Components The momentum equations can be simplified by the introduction of rotary components of horizontal velocity, such that ^ + = 1*,.,+n>,e., 0_ = 1? + I?"; iv , = 0+ - l T , + Ures = U -iv ^ rea re ( D 5 Q ) (D.51) re and the general rotary components E ± (E ,-x±iE ,-y)^ = re re ( D 5 2 ) —* where E Tea is an arbitrary vector, x and y unit vectors in the x and y directions respec- tively. Introducing the surrogate variables (D.50-D.52) into (D.38) removes the Coriolis coupling: —* where the barotropic gradient defined as G (= —g^k'n)- I n the same way, boundary conditions (D.45-D.46) can be written <9i9 ± ^~d7 + T ± A -— + T uz = 1 ± V F m * ± = C W / ^ + r* { z = T , ) ' ( D (z = -h ). m ' 5 4 ) (D.55) Appendix D.7 D. Derivation of the Model 227 Equations Partial Solutions to Momentum Equation The solution to (D.53-D.55) can be written as a combination of two partial solutions {P? and P ± ) : ^(z) = G P (z) + P±{z). ± (D.56) ± 1 —* In a physical sense, the solution is split into a gravity wave part (Pi), where the forcing G is unknown, and a remaining part (P2) where the boundary and body forces are known. Actually, P^ and P are solution to respectively 2 with boundary conditions 8P ± "lh A A l * ( = CWyPi* 1 a v ° = = ' ?7) {z = -h ) m 1 and ,± 0 ( dP A l ± 2 \_ , D ± v , 8T± ± with boundary conditions ,1^ v dz + T± = H^ {Z = TJ), m A l ^ + T* oz = domPt+r* (z = -h ). m Using (D.51-D.52) and (D.56), recovery of ii from i?* is straightforward : 'Pi + p r \ ^ G-ii .(Pi -pr + + 1 1 0 \z x G + (P + P ~)x - i{P+ - P )y. + 2 2 2 (D.59) Appendix D.8 D. Derivation of the Model 228 Equations Momentum Equation Reformulation Equation (D.59) could easily be used to eliminate u* TeB from (D.39). However, the pressure gradient would then appear in a non-standard form. An alternate route is to used the depth-average of (D.56) to eliminate G^. The vertical integration of solution (D.56) gives and therefore iP" — Pi can now be eliminated from (D.56) by using (D.60), such that = - ^ r - - I - -p=t - ~ ' (*)) • J (DM) P 1 Using (D.61), Cioo7. ± _ Q +-0+ &f = o 0 t ± Q ± (D.62) -P , ± where ± _ C r/Pi{z = 1(X 0± = cff? - ^ - P ^ z -h ) (D.63) m (D.64) = -h ). m Hm Using (D.50-D.51), (D.62) can now be written, similarly to (D.59), as ^ = ( ^ ) ^ - i( ^ ) i x ^ - ^ + / J - ) i + ^ - ^ (D-65) Finally, (D.65) can be introduced in (D.39) to eliminate the bottom velocity. The vertically averaged momentum equation then reduces to an equivalent two-dimensional system: —* / ' X U , + Ct'tire, = G , + * : , + R , re re e re +f res " ~77~~) (D.66) Appendix D. Derivation of the Model Equations 229 where all the prime quantities contain contributions from bottom stress (D.67) " *« = ¥, . + ((3 + (3~)x - (D.69) - f3-)y. + re D.9 (D.68) ^r~' = e Generalized Wave Continuity Equation A generalized wave continuity equation (GWCE) can be built by introducing the equivalent two-dimensional form of the vertically averaged momentum equation (D.66) into the vertically integrated continuity equation (D.37). We actually use the so-called weak form of (D.37), derived in Appendix E (E.8) : (<Pi[Vh • (H u .)]) = (w(V • m fc re Q r e t (D.70) )), where ( ) is a domain integral over (x,y) and ipi(x,y) are the basis functions. The left hand side term of (D.70) could be rewritten as From Gauss's theorem, thefirstterm on the right hand side of (D.71) can be expressed as afineintegral over the boundary of the domain, such as, / V • (H u ,<pi)dxdy = / h m re JQ Hu, m re • htfids, (D.72) J8Q where n denotes the unit vector normal to the boundary 50, directed outward and f dxdy expresses the same integral as (). Introducing (D.72) into (D.71), and thereafter Q into (D.70), gives the following weak form of the continuity equation: -(H u . m re • V <Pi) = h ~ I Jd@ Hu m res • fiipids + (<piVh • Qres)- (D.73) Appendix D. Derivation of the Model Equations Introducing (D.66) into (D.73) eliminates u 230 and produces the resulting weak form res of the GWCE : g'gH V <i -fx gH V q (a'Y + f m h res m h r Vh'Pi ) - 2 + a'H A m res -fx / HA m r V <p ) + h (a') + f 2 2 ~ Jd® Hu m ret • rupids (D.74) (<piV -Q ), { h ret where A , includes all the known forcings and nonlinear terms, other than Q , such as re ret — f (D.75) The Galerkin MWR equation (see Appendix E) of the GWCE can be obtained by using the approximate solution JV Sre.{x,y) (D.76) = 53 WfovJi J=l where Cj are the unknown nodal values and N the dimension of the space. Using (D.74) and (D.76), a scalar Helmholtz-like equation in ? , alone is obtained: re [A]{s{z,y) .} re {B(x,y)}-{C(x,y)} a'gH V <pj - f x gH V <pj m h m (a') + / ' 2 a'H K m Bi Ci / Ja® Hu m rea ree -fx • htpids, where [ ] represents a matrix and {} a vector. h 2 Hk m T Vhfi ) + {<fiVfc • Qre.) (D.77) Appendix E Finite Element Method The finite element scheme used to solve the model equations of the thesis (see Appendix D) is based on the Galerkin method of weighted residuals (MWR). An overview of this method follows and details can be found in Celia and Gray (1992). Consider a differential equation that is written using general operator notation as x £ 0, £b{x) = F{x), (E.l) where C is a differential operator, b(x) is the unknown solution of interest, J~(x) is a known forcing function, 0 denotes the domain over which the differential equation applies, and the vector x contains the independent variables. The M W R proceeds byfirstchoosing afinite-dimensionalfunction space, characterized by a set of basis functions {<fi,<f2, • • •, <PN} with AT being the dimension of the space, from which an approximate solution to equation (E.l) is sought. The approximate solution can be expressed as B(x) = Y B < (x), i j Pj (E.2) where Bj is the approximate solution evaluated at the node location Xj (i.e. nodal value) and <pj the basis functions associated with node j. MWR seeks to minimize the amount by which B(x) fails to satisfy the original governing equation (E.l). A measure of this failure can be defined by U(x) = CB(x) - T{x), 231 (E.3) Appendix E. Finite Element 232 Method where TZ(x) is called the residual. This residual is minimized by forcing it to zero in a weighted average sense over the entire domain, i.e. / n(x)1>i(x)dx = 0, Je (E.4) where ibi{x) is a member of a set of weight functions. The Galerkin method chooses weight functions "4>i(x) to be equal to the basis functions <pi(x). Therefore, (E.4) can be written as / (CB - F)<fidx = 0. (E.5) JQ This equation is called the Galerkin MWR equation. The Galerkin method can be seen as a search for a weak solution to the original differential equation. A weak solution of equation (E.l) can be defined as one that satisfies the less stringent equation given by / (Cb)ipi(x)dx = f F{x)<pi{x)dx. JQ (E.6) JQ For example, the weak form of the continuity equation (D.37) can be written / <pi[V • (H u ,)]dxdy = / «#(vX • Q )dxdy, h Je m re ret JQ (E.7) or, for convenience, as (<Pi[Vh • {HjSre.)}) = (>.-(V • <?„..)>, fc (E.8) where ( ) is a symbol used to define the domain contained in the integral. The finite element scheme used in the thesis is composed of 2D (x,y) linear triangular elements and ID (z) linear elements. In the vertical, the basis functions are the piecewise Lagrange polynomials : -- L i tj-zj-l <Pji ) z < Z < Zi - Zi<z<z 0 J i + all other z i ( -) E 9 Appendix E. Finite Element 233 Method Sea S.urface j=3 Z; j=1 J=2 ///////////////////// (a) (b) Figure E.l: Example (a) of the vertical grid and (b) of an horizontal triangular element. where the subscript j refers to node j (see Fig. E.l). In the horizontal, the basis functions are <Pji i y) x = ^ X j + l V j + 2 ~ X J+ ^'+ ) + 2 1 ~ Vi+*) + ( i+2 ~ J+i)v] x x x J ( -!0) E where counterclockwise local node numbering is assumed (see Fig. E.l), indicial additions give 2 + 2 = 3 + 1 = 1 + 3 = 1 and 3 + 2 = 2 + 3 = 2, and A is the area of the triangular element on which the calculation is being done. The area can be expressed as A = \ M 2 / 2 - 1/3) + Z2G/3 - 2/i) + 2:3(2/1 _ 2/2)] • Appendix F Statistical Interpolation This appendix details the statistical interpolation method used in the thesis. A standard reference for this method is Gandin (1965), and more recently Daley (1991). In the oceanographic literature, good references are Bretherton et al. (1976), Denman and Freeland (1985), and Le Traon (1990). F.l Basic Theory Consider the univariate case where the variable to be estimated is solely function of one type of observation. Define i)(r) as the variable, where r indicates the three-dimensional spatial coordinates. Consider ij}(f) as composed of a large-(spatial) scale component (ip(r)) and a mesoscale component (i/'me«o(r))- Now, introduce afirstassumption : Assumption F . l The large-scale component of the variable is known. Using this assumption, an estimate of t/^ri) can be obtained from a linear combination of the N observations : l M ? 0 = Wi) + £ Wiki^in) - Vfe)], where V'ofc* (*"*:) is the observation of the variable at position (F.l) and Wik a weight function indicating the weight given to each observation increment [^obs{rk)—^i^k)] in the analysis. 234 Appendix F. Statistical Interpolation 235 A second assumption is now made : A s s u m p t i o n F.2 The observation values contain errors. where e(fk) is instrumental error and small unresolvable scales of ty. A true value of the variable can then be introduced: i])t{?) = V'(^) + Vme*o(r*)- It is , convenient here to define a simpler notation: if) = ifrettffi), V* — V'C^i)) V'ob* ^obaifk), 7 = eat ipl = ij>t{?i), • • • Now, subtracting from each side of (F.l) gives - ft = ft - ft + E wank. k=i ft. t - ft )- (F.2) 1 The expected analysis error variance (E^) at position r; is found by squaring both sides of (F.2) and applying the expectation operator (Daley, 1991): K = «v>L -ft) )= ((ft -ft) )+ 2E w*< 2 2 + E E WaWuiUk. fc=i J=I ftML - v;)> - ft))- (F.3) Note that the expectation operator commutes with the addition operator. Differentiating (F.3) with respect to each of the weights Wik gives dE 2 ^ N = m L - - ft)(ft - ft)) + 2 E Wu{(i&. - ft){ft<*. - ft))- (F-4) By equalizing equation (F.4) to zero, the weights that minimize E^ axe obtained : E wiMi. - ftML - ft)) = -<(vl - ft)(ft - ft))- (F.5) Appendix F. Statistical Interpolation 236 A third assumption is used in the interpolation method : A s s u m p t i o n F.3 There is no correlation between the observation error and the mesoscale component of the variable. (e(f)Vme.o(f)) = 0 With assumption F.3, we can write that (Ms - ft)(ft - ft)) = -(^Lso) = 0. (F.6) Using (F.6), (Ms - - ft)(ft ft)) = ([(ft^-ft)-(ft-ft)\(ft-ft)) = ((^L - ft)(ft - ft)) - ((ft - tf)(ft - = ft)) -((ft-ft)(ft-ft)) (^meso^meso) (F 7) - and (Ms-ftKftoOs-ft)) = ( [ M s - ^ - i f t - f t M f t ^ s - f t ^ - ^ - f t t ) } ) = (Ms - = (e £ ) + (ft softr so)- k l ft)Ms me - ft)) + (M - ne ft)(ft - ft)) (F.8) The introduction of (F.7) and (F.8) into (F.5) gives E Wu[(ft ft .o) meao me + <eV>] = WLSJLSO) (F.9) or E ^ A i=i W = CW (F.10) Appendix F. Statistical 237 Interpolation (F.10) can also be written as N (F.ll) Wa = i ; C « A i , 1 or N (F.12) It is convenient here to introduce an additional assumption A s s u m p t i o n F.4 The observation and large-scale component of the variable are unbiased. <V>1 - ft) = (e ) = 0 k M - &) = {i> .o) = o me (ft - ft) = (i, J =0 me Using assumption F.4 and the covariance definition of two variables, 00x1(5,17) = (sq) — <-><«> = Cfci = A = cov(ft , w cov(ft , ,ip ), e 0 elo meto V4 , ) + cov(e , e ) = C + cov(e , e ). h l k e 0 l w Using (F.12), (F.l) can now be expressed as N ^stifi) = ftfi) + N £ [=1 CK £ KklMrk) Note that with assumption F.4, ((ft, - ft)) = (ft ba observation increment, ((ip^,, value TJj ,t(ri) e — ft)), - ^(f*)]. (F.13) it=i - ft) - (ft - ft) = 0 and the is thus unbiased. This implies that the estimated defines by (F.13) must also be unbiased. Therefore, (F.13) is an unbiased linear estimate of iff (ri). t The variance of the error in estimating i/'( *i) becomes, inserting (F.5) and (F.7) into T (F.3), K = WL,oftme,o) ~ £ W i k (ft m e ,J m e .o), (F-14) Appendix F. Statistical Interpolation 238 or, using (F.12), 4 = C«-f;C«f;C Ai . 1 K k=l (F.15) 1=1 Cu defines the autocovariance of the mesoscale component of the variable at estimation point Ti. F.2 How to Apply the Theory to a Physical Problem In theory, statistical interpolation is a powerful method but in practice, it could be an hazardous process. Problems with matrix inversion, unknown large-scale component of the variable to be estimated, and unknown covariance structure of the mesoscale component can all be encountered. Firstly, there is the danger that the matrix A_, composed of elements A^, has a very small eigenvalue and then is poorly conditioned. In this case, the matrix inversion would fail. This can happen if two observation stations are very close together. The problem can be avoided by pre-passing through the observations and combining separate observations that are very close together into a single "super-observation", which is thereafter treated as an ordinary observation. Conditioning problems are also substantially improved as the matrix becomes more diagonally dominant. This leads us to the introduction of a fifth assumption : Assumption F.5 Observation errors are uncorrelated with one another. (e e ) — cov(e ,e ) = £8M h l h l where £ is a known variance and Ski the Kronecker delta. The size of the matrix that can be inverted is limited by computer capacity. In many cases, the matrix A_, which has an order equal to the number of observations, is simply too large to invert. The usual solution to this difficulty is to apply the statistical Appendix F. Statistical Interpolation 239 interpolation algorithm locally rather than over the whole domain simultaneously. In other words, each gridpoint is analyzed separately, using only the observations that are likely to be highly correlated with that particular gridpoint. The inversion of one huge matrix is then replaced by the inversion of many smaller matrices. The major difficulty with statistical interpolation lies in the estimation of the necessary ensemble statistics from a severely limited sample of data. In oceanography, ensemble statistics are rarely well-known and experiments are rarely repeated. The large-scale component and the covariance function of the mesoscale component are both usually estimated from a finite sample of realization of observations combined with some a priori prejudices. From Assumption F.l, the large-scalefieldis assumed known. In most of the situations, the best that could be done is to estimate this large-scalefieldfromthe only few available observations. In this case, the mean square error of the estimate is increased by the additional error due to the uncertainty on the large-scale field. For example, a spatially constant meanfieldcan be estimated, from the observations, by (Bretherton et al., 1976) E E »lMrfc) A = ~^n • k E EA/fe ( - ) F 16 1 In this case, (F.15) becomes ( = - 2. ^ 2. ^ fc=1 ' =1 + N N \ ^ 2 - • (F.17) EE » A Note that a spatially variable meanfieldcould also be considered (see Le Traon (1990)). For the covariance function of the mesoscale component, simplifications can be made through additional assumptions: 240 Appendix F. Statistical Interpolation Assumption F.6 The horizontal and vertical structures of the covariance of the mesoscale component of the variable to be estimated are separable. Ckijxkyy^z^x^y^Zj) = C^(a .,yA.,a; -,y )C^(zfc,z ) A t t 1 Assumption F.7 The horizontal covariance of the mesoscale component of the variable is homogeneous. C^(xfc, yjfe, xi, (xi — Xk,yi — yk) yjfe) depends only on the relative horizontal displacement rather than on the absolute horizontal locations = » ki( k,yk,Xi,y ) C x k = C^(z,- - x ,yj k (xi,yi) and (xk,Vk) - y ). k The problem is now to choose a suitable covariance model for the mesoscale component. Two constraints are applied to this choice (Denman and Freeland, 1985) : • ^S"** A N D 9 3 C 8S~ f i ) are continuous at \f k - TJ| = 0. • The spectral densities of C(fk — r») are strictly positive. This constraint implies that C(fk — Ti) —> 0 when |f* — fj| 3> Ld, where Ld is some natural length scale of thefieldbeing sampled.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modelling of the circulation of northern British Columbia...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modelling of the circulation of northern British Columbia waters Jacques, Renee 1997
pdf
Page Metadata
Item Metadata
Title | Modelling of the circulation of northern British Columbia waters |
Creator |
Jacques, Renee |
Date Issued | 1997 |
Description | A numerical model is used to improve our understanding of the mean circulation of northern British Columbia waters. The velocity field arising from tidal, wind, river runoff, and baroclinic effects is studied and the influence of the earth's rotation investigated. The model is three-dimensional, diagnostic, nonlinear through tidal forcing, and uses finite element techniques. The grid has linear horizontal triangular elements with vertical sigma levels. The model output is a linear combination of the currents obtained under the individual forcings. In this thesis, the model has been modified in order to improve the parameterization of friction processes, to add the vertical advection of tidal momentum by the tides, and to include the stratification effect on eddy viscosity. A major circulation feature is a large cyclonic eddy covering the central-eastern Dixon Entrance. This eddy, referred as the Rose Spit Eddy, is actually simulated by a combination of two or three smaller gyres, which seasonally vary in size and intensity. The modeled eddy is, in good agreement with observations, characterized by a rotation period of 3 to 40 days, depending on the season and the initial distance of the numerical particle from the center of the gyres. This eddy persists in all seasons but occasionally disappears. A disappearance occurred, consistent with current observations, in October 1984. The Learmonth Bank Gyre, an anticyclonic eddy located in western Dixon Entrance, is also simulated by the model. Baroclinicity and tidal rectification, respectively responsible for the generation of the Rose Spit Eddy and Learmonth Bank Gyre, are the dominant effects in the study region. The earth's rotation, necessary for the formation of those two eddying motions, as well as the vertical advection of tidal momentum were found to be essential to bring the current speeds closer to the observed values. The effect of stratification on the residual currents, through turbulence reduction, was found to be negligible. |
Extent | 11567471 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0053281 |
URI | http://hdl.handle.net/2429/6696 |
Degree |
Doctor of Philosophy - PhD |
Program |
Oceanography |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-250725.pdf [ 11.03MB ]
- Metadata
- JSON: 831-1.0053281.json
- JSON-LD: 831-1.0053281-ld.json
- RDF/XML (Pretty): 831-1.0053281-rdf.xml
- RDF/JSON: 831-1.0053281-rdf.json
- Turtle: 831-1.0053281-turtle.txt
- N-Triples: 831-1.0053281-rdf-ntriples.txt
- Original Record: 831-1.0053281-source.json
- Full Text
- 831-1.0053281-fulltext.txt
- Citation
- 831-1.0053281.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0053281/manifest