## Abstract

The tidal response of northern British Columbia coastal waters is studied through simulations with a three-dimensional, prognostic, primitive equation model. The model is forced at the boundaries with the leading semidiurnal and diurnal constituents and experiments with stratified and homogeneous fluid are compared. The barotropic response shows good agreement with previously published studies of tides in the region. A comparison with tide gauge measurements indicates that average relative rms differences between observations and the model surface elevation field are less than 5% for the largest constituents.

An internal tide is generated in cases where the model is initialized with a vertical stratification. Diagnostic calculations of the baroclinic energy flux are used to identify regions of generation and propagation of internal tidal energy. With a representative summer stratification, the integrated offshore flux is about 0.5 gigawatts, higher than previously estimated from theoretical models. Comparisons between observed and modeled *M*_{2} current ellipses are discussed for several moorings and demonstrate the significant influence of the internal tide.

## 1. Introduction

Motions at tidal frequency in the coastal ocean are often composed of an astronomically driven surface (barotropic) tide and an internal (baroclinic) tide, arising from interaction of the barotropically forced, stratified fluid with bathymetry. In contrast to the barotropic tide, the baroclinic component contributes relatively little to oscillations of the sea surface; its signature appears primarily in tidal currents and internal isopycnal motions. Numerical models have been applied in several instances to study the response of a coastal region to specified tidal forcing. To date, however, models of open coastlines have been concerned almost exclusively with the barotropic tide. (See Foreman et al. 1995 for a review.) While baroclinic models in three dimensions have been applied to embayments and channels of limited extent (Matsuyama 1985; Wang 1989), to our knowledge, baroclinic tidal studies of open coastlines have been limited to applications of two-dimensional, cross-shelf models (e.g., Craig 1988; Sherwin and Taylor 1990). Recently, Holloway (1996) applied an essentially two-dimensional version of the Princeton Ocean Model to the Australian North West Shelf. Results, particularly *M*_{2} baroclinic energy fluxes, were highly sensitive to specification of the cross-shelf topographic section.

Tides along the northern coast of British Columbia have been the subject of a number of previous modeling studies. Flather (1987, hereinafter F87) forced a finite difference barotropic model of the northeast Pacific with the *M*_{2} and *K*_{1} tidal constituents. Foreman et al. (1993, hereinafter FHWB93) presented results from a barotropic finite element model with eight constituents for a somewhat smaller region. Extensive comparisons with sea level data in these two papers showed generally excellent agreement between observations and the model surface elevation field. This work was extended recently by Ballantyne et al. (1996, hereinafter BFCJ96) who employ a three-dimensional diagnostic finite element model. The performance of the three-dimensional model with respect to sea level is comparable to the barotropic models. In addition, detailed comparisons of tidal current ellipses with observations through the water column were presented for various locations. For the *M*_{2} constituent, significant discrepancies were found that may be attributable to the presence of baroclinic tides. In each of the aforementioned studies, the nature of the model precluded the possibility of an internal tidal response.

In the present paper, tidal simulations of the northern coastal waters of British Columbia with a fully prognostic three-dimensional model are discussed. The numerical formulation is based on the version of the Princeton Ocean Model (G. Mellor 1993, personal communication) used by Oey and Chen (1992) to simulate flows over the northeast Atlantic shelves. The model is driven at the boundaries with the principal diurnal and semidiurnal constituents. Numerical experiments with homogeneous and with vertically stratified fluid are compared. With homogeneous fluid, the solution is similar to that of BFCJ96. The barotropic tide is well represented, but current ellipses in many locations bear little resemblance to observations. Inclusion of a vertical stratification allows an internal tide to be generated, significantly affecting current ellipses for the semidiurnal constituents. This leads to an improved representation of *M*_{2} currents in several regions over the shelf and slope. Calculations of the baroclinic energy flux are used to identify regions of internal tide generation and propagation.

In the next section, a brief review of the regional oceanography and geography is provided. This is followed by a section containing a description of the model configuration and outlining the numerical experiments. Section 4 discusses the barotropic tidal response of the model, while the baroclinic response is described in section 5. A summary and suggestions for future work are given in the concluding section.

## 2. Regional oceanography

A thorough review of the physical oceanography of northern coast of British Columbia is available in Thomson (1989). Our purpose here is to draw attention to the features of the region that are pertinent to discussion of the numerical simulations. The coastal geography and bathymetry of the region are illustrated in Figs. 1a and 1b. The bathymetry is characterized by a broad, *O*(100 km wide) shelf, which is partially sheltered from the open ocean by the presence of Moresby and Graham Islands, known collectively as the Queen Charlotte Islands. The shelf is narrow on the west side of these islands with depths dropping to 2500 m over a distance no greater than 30 km from the coast. Three connected seas separate the Queen Charlottes from the coastal mainland. Queen Charlotte Sound is open to the Pacific on its western side and marked by three relatively shallow banks. Proceeding from south to north, these are Cook Bank, Goose Island Bank and Middle Bank. To the north, Hecate Strait is a broad, shallow passage with an elongated submarine valley on its eastern side. Owing to the presence of Dogfish Banks, Hecate Strait is very shallow (less than 30-m depth) on its western side. The third significant sea is Dixon Entrance, a passage of O (60 km) width connecting northern Hecate Strait to the open Pacific. Learmonth Bank is a shallow barrier found near the opening of Dixon Entrance.

Barotropic tides of the semidiurnal and diurnal frequency bands propagate approximately as large-scale Kelvin waves in a cyclonic sense around amphidromes in the North Pacific. Amplitudes of both bands are significant in the northeast Pacific, and hence tides along the coast of British Columbia are classified as mixed, predominantly semidiurnal. Mean tidal ranges are about 2 m in deep water, increasing to about 3.5 m near the entrances of Queen Charlotte Sound and Dixon Entrance and to about 5 m near Prince Rupert. The tides are subject to a spring–neap cycle, which arises through an interference of *M*_{2} and *S*_{2}, the dominant semidiurnal constituents, and also *K*_{1} and *O*_{1}, the largest diurnal constituents.

In the presence of a sloping bottom, barotropic tides can force vertical oscillations of a stratifed fluid and thus provide a source of baroclinic energy. In two dimensions (*x, z*), the dynamics can be expressed by an equation governing the streamfunction (Huthnance 1989) *ψ*,

where *Q* is the barotropic cross-shelf volume flux, *H* the bottom depth, and *ω* the tidal frequency; *N* = (−*g**ρ*_{z}/*ρ*_{0})^{1/2} is the buoyancy frequency with *ρ* the density, *ρ*_{0} a reference density, and *g* the gravitational constant. Provided that (*ω*^{2} > *f*^{2}), (1) is hyperbolic with propagation along rays whose slope is given by

where *f* is the local Coriolis parameter. The bottom slope *α* is supercritical if *α*/tan*c* > 1, and propagation toward deep water is anticipated. Conversely, a subcritical slope is associated with onshore propagation.

Using typical temperature and salinity data from the continental slope region to determine *N* (section 3, Fig. 3b), the ratio *α*/tan*c* is found to easily exceed unity over the slope. Stratification in the upper 100–200 m of the water column undergoes a seasonal modulation, associated with coastal runoff and surface warming during summer months and vertical mixing during winter. Below 200 m such variations are weak and internal tide generation with offshore propagation thus appears possible throughout the year over the slope. Although the barotropic forcing is constant, the internal tide may nevertheless be intermittent as physical processes affecting the stratification (e.g., mixing or wind-induced downwelling) can modulate generation or propagation. Drakopoulos and Marsden (1993) examined evidence for internal tides off the west coast of Vancouver Island, to the southeast of our study area. They observed offshore propagation and a seasonal modulation, with maximum intensity during summer months.

## 3. Numerical model

The numerical model used in this study is based on the Princton Ocean Model (G. Mellor 1993, personal communication), and a detailed description is given in Oey and Chen (1992). Briefly, the model solves finite difference analogs of the primitive equations, with a nonlinear equation of state relating density to the temperature and salinity. The level-2.5 turbulent kinetic energy closure scheme of Mellor and Yamada (1982) is used to determine vertical mixing coefficients for momentum and tracer (temperature and salinity) fields. Horizontal mixing coefficients are modeled with the Smagorinsky formulation for diffusion. Bottom frictional stresses are modeled by a quadratic drag law with a drag coefficient formulated to yield a logarithmic bottom boundary layer. Buoyancy and momentum fluxes through the surface have been set to zero.

The numerical algorithm is described by Mellor in the user’s guide (1993) for the Princeton Ocean Model. The model employs a C grid in the horizontal and a terrain following sigma coordinate system in the vertical. A mode-splitting technique is applied for integration of external and internal modes with different time steps. The integration proceeds according to an explicit leapfrog time stepping scheme, except for terms involving vertical diffusion of momentum and that are treated implicitly. A weak time filter is included to prevent development of a computational mode with the leapfrog scheme.

Through a forced gravity wave radiation condition (see Oey and Chen 1992; F87), the model is driven at open boundaries with the four leading tidal constituents. Boundary elevation data for these constituents are obtained from the global tidal model of Schwiderski (1979, 1981a–c). Depth-averaged tidal velocities required at open boundaries for the radiation condition are obtained from an iterative procedure similar to one outlined in F87. Boundary conditions differ from those of Oey and Chen (1992) in certain respects. An Orlanski radiation condition has been applied to the internal mode velocities, temperature, and salinity at the open boundaries.

Barotropic tidal motions are driven entirely by oscillations at open boundaries; the earth tide, load tide, and tide-generating potential have been neglected. With a coastal model of similar extent, FHWB93 concluded that these forcings produce only slight modifications to semidiurnal constituents.

### a. Model grid

The governing equations of the model are solved over a Cartesian grid with a uniform resolution of 5 km in the horizontal (*x, y*) coordinates. The universal transverse mercator map projection, used by Hannah et al. (1991), relates the (*x, y*) coordinate of each grid point to latitude–longitude coordinates. A rotation of 30° is applied to align the *y* coordinate of the model approximately with the shelf break, thus reducing the number of land points. The Coriolis parameter is set to a constant appropriate to 53°N. In most cases, there are 21 sigma levels in the vertical with the following distribution: *σ*_{k} = (0, −0.003, −0.006, −0.015, −0.03, −0.05, −0.07, −0.09, −0.11, −0.13, −0.165, −0.2, −0.3, −0.4, −0.5, −0.6, −0.7, −0.8, −0.9, −0.95, −1), where *σ*_{k} = 0 at the surface and *σ*_{k} = −1 at the bottom. Velocity and tracer fields are defined at midpoint between these levels.

The model topography was constructed by averaging depth data within a 5 km radius at each grid point. The depth data for most of the domain were obtained from a high resolution (1 km) topographic database available for the British Columbia coast. It was necessary to supplement this with ETOPO5 topography for the northern reaches of Dixon Entrance and for regions located well offshore of the continental slope. No smoothing was done to the resulting topographic field except for an area adjacent to the offshore boundary (0 < *X* < 150 km), where a smoother was applied, eliminating the presence of Bowie and Union seamounts from the model. The resulting model domain is shown in Fig. 2 with contours of bottom topography. The domain is open on the two cross-shelf boundaries and on the offshore boundary. Two narrow straits, Johnstone Strait to the south and Clarence Strait to the north (Fig. 1) have been closed off.

### b. Experiments

The numerical experiments of this study are listed in Table 1. Most of the discussion is focused on a comparison of the first two cases, but mention is also made of some additional runs. The first case, denoted 14v, is a control run in which a homogeneous unstratified fluid is specified. The baroclinic response of the model was examined with a second experiment, denoted 14u. Here the temperature and salinity are initialized as horizontally homogeneous vertically stratified fields representative of summer stratification. Figure 3a illustrates the vertical variation of temperature and salinity used for experiment 14u. The vertical profile is specified with an analytic form composed of a linear region near the surface and an exponential form below. This profile, although idealized, provides a reasonable approximation to the density stratification typically observed over the continental slope and outer shelf during summer conditions (Fig. 3b). Additional experiments include a case with a “winter” stratification and a stratified diagnostic case. The initial stratification for the winter case, 14s, differs from the profile of Fig. 3a only in the upper 125 m where, to allow for a winter mixed layer, the fluid is homogeneous.

Motivation for the diagnostic (fixed temperature and salinity) experiment, 14p, is to examine sensitivity to stratification in the absence of an internal tide. It is known that properties of tidal current ellipses can be sensitive to the value of the vertical viscosity (Prandle 1982). In the present case, this parameter is determined as part of the numerical solution through the turbulence closure submodel. Its value will depend on the several factors, in particular, the vertical stratification. In experiment 14p the fluid was stratified as in 14u, but an internal tidal response was precluded by relaxing the temperature and salinity fields to their initial values with a very short time constant. Except in the bottom boundary layer where relatively minor differences were obtained, results from this diagnostic case, in particular tidal ellipse parameters, were nearly indistinguishable from those of 14v. The results of experiment 14p demonstrate that the significant differences between 14u and 14v discussed below are not merely a result of stratification-dependent viscosity, but rather the presence of internal tides in 14u.

For experiments 14v and 14u, the model was integrated for 34 days of simulated time. Hourly fields from the last 29 days were subject to an harmonic analysis to separate the response at the four constituent frequencies. This analysis period is sufficient for the harmonic analysis to separate the response within the diurnal and semidiurnal frequency bands, but also short enough that the mean stratification does not change appreciably. The initial spinup time is sufficient to establish an equilibrium for the tide.

### c. Resolution tests

Simulations of the baroclinic tide require resolution of the scales of the topography and the wavelength associated with the internal tide. It is also desirable to minimize pressure gradient error and other numerical truncation errors associated with sigma levels (Mellor et al. 1994). A length scale, *H*/tan*c,* based on (1) (Huthnance 1989) yields a value of 25 km over the shelf, assuming *H* = 200 m and *N* = 1 × 10^{−2} rad s^{−1}. In 2000 m of water with a typical deep value for *N* (3 × 10^{−3} rad s^{−1}), the length scale is 75 km. The 5-km grid resolution appears sufficient for deep regions, but may be marginal over the shallower regions of the shelf.

A limited convergence study was done to examine sensitivity of the baroclinic response to variations in horizontal and vertical resolution. These cases are the last three entries in Table 1 and include only forcing with the *M*_{2} constituent. The runs were restricted to a duration of 5 or 6 days with harmonic analysis applied to the last 3 or 4 days. In experiment RT1 the basic 5-km grid with 21 vertical levels is maintained. The vertical resolution is refined in experiment RT2 to 41 levels, while in RT3 the horizontal resolution is doubled to 2.5 km. The additional vertical levels in RT2 were placed at midpoint between those of the 21-level distribution. The topography in RT3 was obtained from linear interpolation of the 5-km grid and does not introduce finer topographic scales.

An intercomparison of results from these three cases showed only minor quantitative differences in the baroclinic currents of the model to variations in the horizontal or vertical resolution. The radiation pattern of the internal tide was essentially unchanged in these cases and resembled the *M*_{2} response of the stratified case with four constituents (14u). Thus, the model does not appear to be very sensitive to changes in resolution. It is likely that a model that resolved finer topographic scales would introduce new finestructure to the internal tidal field. However, as shown below, the large-scale structure of the internal tide is determined by the larger scales of the topography and the barotropic tide, which are both well resolved.

## 4. Barotropic response

The results in this section are drawn mainly from experiment 14v with homogeneous fluid. Stratification has little effect on the barotropic response of the model. There are, however, subtle differences between 14v and the stratified runs, the most notable of which is an additional loss of barotropic tidal energy over the shelf with stratified fluid.

### a. Surface elevation

The amplitude and phase of the surface elevation field at the *M*_{2} and *K*_{1} frequencies are presented in Figs. 4a and 4b. The *M*_{2} phase progresses steadily northward along the outer coast. As the tide propagates into Hecate Strait from the south, amplitudes increase due to shoaling effects on the shelf and reach almost 2 m on the east side of Graham Island. The shoaling also induces an east–west phase difference of about 10° across the Strait. The *M*_{2} amplitude and phase of Fig. 4a agree closely with comparable fields presented in earlier studies (F87, FHWB93, BFCJ96). The most notable difference with these last two studies is a small discrepancy in the phase over the northern part of the domain. In Fig. 4a, the 270° phase contour abuts on the outer coast of Alaska, whereas it passes through northern Dixon Entrance in the finite element models. One factor that may account for this difference is that FHWB93 and BFCJ96 tuned their boundary data (also taken from the Schwiderski atlases) to obtain an improved agreement with sea level data over their domain.

Figure 4b shows that the sea level amplitude field for *K*_{1} is about 40% of *M*_{2} in the deep ocean with relatively little amplification in shallow regions. The amplitude and phase distribution for *K*_{1} also agrees with earlier studies, even with respect to small-scale maxima and minima over the shelf. Some discrepancies with respect to the distribution of the phase are evident over the deep ocean regions where our results are in closer agreement with those of F87 than either FHWB93 or BFCJ96. Once again, a likely reason for this is the tuning of boundary data in these studies. (F87 made no such adjustments to his boundary forcing.)

Distributions for the S_{2} and *O*_{1} constituents (not shown) are similar in pattern to *M*_{2} and *K*_{1}, respectively, and generally agree with unpublished results, provided by M. Foreman, from the finite element models. Here S_{2} amplitudes are about 30% of *M*_{2}, while *O*_{1} amplitudes are about 60% of *K*_{1}.

To provide a quantitative assessment of the elevation field, averages of the absolute rms error, *L*^{−1} Σ_{L }*D,* and the relative rms error, *L*^{−1} Σ_{L }*D*/*A*_{o}, were computed for each constituent. The averaging is based on observations from 50 (=*L*) tide gauge sites, whose distribution is shown on Fig. 1a; *D* is the rms difference over a tidal cycle between model and observations, given by

where *A* and *ϕ* are amplitudes and phases for a given constituent, and subscripts *m* and *o* refer to model and observations, respectively. The 50 gauge sites are a subset of the 63 listed in Table 2 of FHWB93. As they are mainly distributed around the shelf region, where the greatest variations occur, these sites provide adequate coverage of elevation field. For a given tide gauge, the closest corresponding model grid point was chosen for comparison. Of the 13 gauges omitted from consideration, 12 were situated at locations that were unresolved by the 5-km grid, usually because they were at the head of narrow inlets. One site (Alert Bay) was omitted because it was close to the closed-off entrance to Johnstone Strait (Fig. 1a), where realistic results are neither expected nor obtained.

The average rms elevation errors from the comparison are given in Table 2 for the two basic cases, 14v and 14u. The average relative error is less than 5% for each constituent except *O*_{1}, which is slightly higher. As expected, the two experiments have similar differences with observations. The largest absolute errors (e.g., *D* ≈ 24 cm for *M*_{2}) occur near the artificially closed entrance to Johnson Strait. However, the region of large error is locally confined and errors for the nearest available stations around Queen Charlotte Strait are generally smaller than the mean values of Table 2. Based on these results, we conclude that the barotropic tide is reasonably well represented in the model.

### b. Tidal current ellipses

Tidal current ellipses for the depth-averaged velocity are presented in Fig. 5 for the *M*_{2} constituent. In deep regions, tidal currents are small, *O*(2 cm s^{−1}), and it is only over the shelf that the ellipses assume a significant magnitude. Large clockwise-rotating ellipses with relatively large eccentricity (ratio of semiminor to semimajor axis) are found over the three shallow banks in Queen Charlotte Sound. Particularly large ellipses also occur in the vicinity of Cape St. James. In Hecate Strait and especially in Dixon Entrance, the channel boundaries constrain the motions to be nearly rectilinear, aligned along the axes of these channels. Near Rose Spit, the transition point between the two channels, there is a rapid rotation in the orientation of the ellipses. FHWB93 show *M*_{2} tidal ellipses for northern Hecate Strait, which are in good agreement with Fig. 5. In comparison to *M*_{2}, ellipses for *K*_{1} (not shown) are notably weaker in amplitude, especially in Dixon Entrance. The largest *K*_{1} currents are found off Cape St. James, over the shallow banks in Queen Charlotte Sound and on the Alaskan shelf.

With homogeneous water, the tidal ellipses show little depth dependence above the bottom Ekman layer. Figure 6 compares the vertical variation of *M*_{2} and *K*_{1} tidal ellipses from experiment 14v with observations at station W02 located in central Hecate Strait (see Fig. 1a). This comparison, and others from Hecate Strait, show that a homogeneous motel captures the essential aspects of *M*_{2} and *K*_{1} tidal currents in this region. Except in the bottom boundary layer, tidal ellipses at W02 in the stratified experiments (not shown) are virtually identical to those of Fig. 6, suggesting that internal tidal motions are absent at this site. BFCJ96 compared model and observed ellipses for four sites, including W02. At each of these sites, their results are in close agreement with tidal ellipses computed from our experiment with homogeneous fluid. However, it is only at the W02 site that good agreement with observations is obtained. At the other sites, located in Queen Charlotte Sound, Dixon Entrance, and on the west side of the Queen Charlotte Islands, observations show significant vertical structure that the model of BFCJ96 and our homogeneous fluid model cannot capture.

### c. Energy fluxes

Propagation of the barotropic tide may be examined from consideration of depth-integrated energy fluxes associated with each tidal constituent. These fluxes are determined from the harmonic analysis of the barotropic currents and the surface elevation field according to a method outlined in appendix A. The resulting energy flux vectors are shown in Fig. 7 for the principal semidiurnal and diurnal constituents. In both cases, the field is dominated in the deep water by the northward propagation of energy associated with large-scale gyres of tidal energy in the Pacific. Relatively little *K*_{1} energy enters Hecate Strait or Dixon Entrance. Near the shelf break, the *K*_{1} constituent has circular gyres of energy flux located south of Cape St. James and along the Alaskan shelf. F87 suggested that similar features appearing at these locations in his solution were associated with a diurnal shelf wave component.

For the *M*_{2} constituent, barotropic energy leaks onto the shelf in a region of O (40 km) width south of Cape St. James. A fraction of this energy continues northward through Hecate Strait and Dixon Entrance, eventually rejoining the deep water flux. A second branch recirculates southward over the shelf, crossing over Goose Island Bank and Cook Bank, where it then rejoins the deep water flux. There are significant *M*_{2} energy fluxes across the shelf break at three distinct locations: along the northern shelf region off Dixon Entrance, south of Cape St. James, and farther south off Cook Bank. A gap, where relatively little energy flux crosses the shelf break, exists between the latter two locations.

The *M*_{2} energy flux was integrated across various transects to infer the rate of energy loss in the barotropic tide over the shelf. The mean integrated flux of barotropic energy into an enclosed region gives the net rate of energy loss to all dissipation mechanisms, of which bottom friction is the most significant. In stratified experiments, conversion to baroclinic energy represents an additional potential sink for barotropic energy. Stratification effects may also modify bottom currents and hence the dissipation rate associated with bottom friction.

Figure 8 illustrates three subregions of the model and gives the net fluxes across the bounding transects for experiments 14v and 14u. Barotropic losses within a subregion are the residual obtained from summing fluxes across the bounding transects. These values are given in Table 3 for the two experiments. In the unstratified case, 14v, a net flux of about 4 GW onto the southern shelf from the deep water is largely balanced by the flux of 3.5 GW into Hecate Strait (1 GW = 10^{9} watts). Hence there is a mean dissipation rate of about 0.5 GW over the Queen Charlotte Sound shelf. Most of the dissipation, about 1.2 GW, occurs in Hecate Strait/Dixon Entrance where relatively shallow water is found.

In the stratified case (14u), the northward flux into Hecate Strait and out of Dixon Entrance is similar to the homogeneous case and the rate of energy loss in this region is essentially unchanged. However, the net flux onto the Queen Charlotte Sound increases to 4.9 GW, implying an additional loss of about 1 GW over the southern shelf. Overall, losses in the *M*_{2} tide over the shelf increase from 1.85 to 3.05 GW with introduction of the stratification. As discussed in the next section, baroclinic flux calculations suggest that conversion to internal tidal energy accounts for a significant fraction of the additional barotropic losses.

## 5. Baroclinic response

To examine the baroclinic response of the model, we consider results primarily from experiment 14u, the basic stratified case. As indicated above, the barotropic response is nearly (but not completely) unchanged in this run from the unstratified case, 14v. However, a vigorous internal tide is generated in 14u that affects tidal currents over much of the domain.

The regions of internal tidal influence are delineated in Fig. 9, which shows contours of the rms *M*_{2} baroclinic current on a near-surface velocity sigma level. Here the baroclinic current is defined as the difference between the total and the depth-averaged current, as in Eq. (A5). The intensity of the near-surface baroclinic component of the tidal current gives an indication of the strength of the internal tide. (The comparable field from 14v is close to zero everywhere; with homogeneous water, near-surface and depth-averaged currents are nearly identical.) The results of Fig. 9 indicate that the internal tidal influence is greatest over the outer shelf of Queen Charlotte Sound and on either side of Learmonth Bank in Dixon Entrance. There is also a significant offshore influence emanating from the shelf break of these two regions. In contrast, the internal tidal signal is relatively weak for the inner portions of the shelf and in Hecate Strait. In addition, a region of comparatively weak internal tides is indicated in deep water off the west coast of the Queen Charlotte Islands.

In the deep water (*H* > 2500 m) beyond the shelf, rms *M*_{2} currents from experiment 14u near the surface (*σ* = −0.0105) have an average value of 6.5 cm s^{−1} with a standard deviation of 3.3 cm s^{−1}. This compares with an average of 2.1 cm s^{−1} (standard deviation 0.3 cm s^{−1}) for the rms depth-averaged *M*_{2} current. Near-surface currents over deep regions are thus generally dominated by the baroclinic component. These results are consistent with values that Thomson et al. (1997) obtained from an analysis of near-surface drifters deployed in the region of Ocean Weather Station P (50°N, 145°W) and propagating eastward toward British Columbia. Their average rms semidiurnal currents were practically the same for both 15-m and 120-m drogued drifters (5.9 cm s^{−1}, with a standard deviation of 2.1 cm s^{−1}). The minimum observed rms semidiurnal values (2.9 cm s^{−1}) likely correspond to the barotropic currents.

### a. Internal motions

The baroclinic response of the model is characterized by internal vertical displacements and baroclinic currents at the constituent frequencies of the barotropic tide. The internal displacements can be examined by decomposing the vertical velocity into a series of harmonics of the form *w*_{n}cos(*ω*_{n}*t* − *λ*_{n}) where *w*_{n} and *λ*_{n} are the amplitude and phase at frequency *ω*_{n}. Since *w* = *d**ζ*/*dt,* where *ζ* is the vertical displacement, the amplitude and phase of the latter is given by *w*_{n}/*ω*_{n} and (*λ*_{n} + 90°), respectively.

Figure 10a presents contours of the vertical displacement amplitude at *M*_{2} frequency for a representative *x*–*z* cross section that intersects the continental margin at Dixon Entrance. Semidiurnal *M*_{2} amplitudes in excess of 40 m are obtained over the steepest parts of the continental slope. These large values result from the direct forcing of vertical motions over the slope by the barotropic tide and indicate generation regions for the internal tide. Here the topographic slope is steeper than the slope of characteristics (*α*/tan*c* > 1) and propagation is directed toward deep water. In the deep regions, displacement amplitudes in the upper 500 m range from 1 to 3 m near the surface to a maximum of about 10 m at mid depth. These near-surface values are comparable to semidiurnal displacements of 1–4 observed in the upper 200 m at station P (Thomson and Tabata 1982). Displacement amplitudes for *S*_{2} (not shown) are similar in pattern to the *M*_{2} field but smaller by a factor of ∼3 along this section.

Contours of vertical displacement amplitude along the cross section for *K*_{1} are presented in Fig. 10b. At the *K*_{1} frequency, significant (>15 m) displacements are manifest over the continental slope. However, at the latitude of the B.C. coast, diurnal frequencies are smaller than the inertial frequency (*ω* < *f*). In this frequency range Eq. (1) ceases to be hyperbolic, and free internal wave propagation is no longer possible. As a result *K*_{1} motions are trapped over the slope and there is no radiation of energy offshore. A similar trapping occurs with the *O*_{1} field (not shown).

A cross-sectional view of the rms *M*_{2} baroclinic current is given in Fig. 11. In deep water the vertical structure resembles that of a baroclinic internal wave mode. There are maxima near the surface with rms values reaching 5–10 cm s^{−1}, minima at approximately middepth, and maxima of 1–2 cm s^{−1} at abyssal depths.

### b. Energy fluxes

Baroclinic energy fluxes can provide insight into generation and propagation of the internal tide. Three different views of this energy flux are presented in Figs. 12a–c. Figure 12a shows, for the entire model domain, the depth-integrated, baroclinic energy flux vectors at the *M*_{2} frequency. (See appendix A for details.) Close-up views for northern and southern subregions of the shelf are included as Figs. 12b and 12c. Energy fluxes for *S*_{2} (not shown) are similar in pattern to those for *M*_{2} but with amplitudes reduced by a factor of ∼6. The baroclinic fluxes are generally much smaller than the corresponding barotropic fluxes. For example, the *M*_{2} energy flux represented by the scale vector in Fig. 12a is 1% of that in Fig. 7a.

The vector field of Fig. 12a is dominated by an offshore flux of baroclinic energy originating from the slope region near the 1000-m contour. The field is highly inhomogeneous with offshore fluxes generated at three principal locations. Along the shelf edge between Cape St. James and Vancouver Island, two generation regions are obtained, separated by a length of shelf where the offshore flux is comparatively weak. The southernmost of these two source regions lies off Cook Bank. Farther north, there is a stretch of the shelf break 150 km long off Dixon Entrance and Alaska where a large offshore flux of baroclinic energy is generated. Also evident is a region off the west side of the Queen Charlotte Islands where the offshore flux is relatively weak. It is this “shadow zone” that accounts for the relatively weak near-surface signal obtained off the Queen Charlottes in Fig. 9. In the far field, offshore fluxes spread geometrically from source regions at the slope.

Generation sites for baroclinic energy in the deep water correspond precisely to the three regions mentioned in section 4c where significant fluxes of barotropic energy cross the shelf break. In these regions the cross-slope component of the barotropic *M*_{2} tidal ellipses (Fig. 7a) is largest and stratified fluid is forced most vigorously up and down the continental slope. In terms of the two-dimensional model [Eq. (1)], these are regions where the cross-slope volume flux *Q* is maximum. On the west side of the Queen Charlottes a different situation prevails; cross-slope motions there are weak due to narrowness of the shelf and presence of the land mass.

An inhomogeneous pattern of baroclinic *M*_{2} energy generation and propagation over the shelf is evident in Figs. 12b and 12c. South of Cape St. James, onshore fluxes into Queen Charlotte Sound are obtained, which originate irregularly from near the 300-m depth contour over a distance extending for about 100 km (Fig. 12b). Near Cape St. James, some of this baroclinic energy leaks into the southwestern reaches of Hecate Strait. There is also a broad flux of baroclinic energy radiating off the southern flank of Goose Island Bank, near the 100-m contour, toward northern Vancouver Island. The pattern of energy fluxes in the region of Dixon Entrance (Fig. 12c) is particularly inhomogeneous. It is possible to discern both onshore and offshore fluxes off the sides of Learmonth Bank. The results also show fluxes directed in a cross-strait sense from the shoaling topography near the northern Queen Charlotte Islands.

Since a large fraction of the baroclinic energy is direct offshore, an estimate of the net baroclinic energy generation may be obtained by integrating the flux normal to transects forming a box around the shelf region. The result of such a calculation for the *M*_{2} and S_{2} constituents is given in Table 4 for different stratified experiments. The integrated baroclinic flux is reduce somewhat from 14u in the winter stratification case (14s), but the basic pattern is essentially unchanged. Offshore fluxes in the three resolution test cases (RT1, RT2, and RT3) show mutual agreement and differ only slightly from 14u, confirming that energy fluxes are relatively insensitive to variations in horizontal or vertical resolution.

Of the net offshore baroclinic flux of 0.56 GW in the basic case (14u), 0.17 GW emanates from the Dixon Entrance/Alaskan shelf region, while 0.34 GW originates from Queen Charlotte Sound. The values represent the baroclinic energy radiating offshore but do not include fluxes radiating shoreward from the shelf break or local fluxes over the shelf or slope. It was noted in section 4 that barotropic losses over the northern shelf region were augmented by 0.18 GW and over the southern shelf region by 1.03 GW. Hence, baroclinic energy conversion readily accounts for the additional barotropic losses over the northern shelf region and a significant fraction of the additional losses over Queen Charlotte Sound. Bottom friction over shelves is the main dissipation mechanism for the barotropic tide and additional losses with stratification may be partly attributable to enhancement of bottom dissipation. Rms currents above the bottom are significantly increased due to the internal tide in some regions over shelf in Queen Charlotte Sound and over the adjoining slope.

The offshore *M*_{2} baroclinic energy flux of about half a gigawatt obtained here is significantly larger than an earlier estimate made by Baines (1982) with a two-dimensional analytical model. Values for the northeast Pacific were not explicitly given but were implied to be less than 0.1 GW. Reasons for the discrepancy are not clear. However, our results show localized sources and the inhomogeneous nature of baroclinic energy generation, features that may be difficult to account for with a two-dimensional model.

### c. Comparison of *M*_{2} current ellipses with observations

As a means of verifying the model, comparisons between observed *M*_{2} tidal ellipses and results from experiments 14v and 14u are presented. Current meter moorings from regions where an internal tidal component is indicated are considered. The discussion is limited to *M*_{2}, the leading semidiurnal constituent, however, ellipses for *S*_{2} are similarly affected by the stratification. Except for station QE6, the observations are drawn from summer deployments (nominally mid-May to Sep) made during the late 1970s and early 1980s as part of the CODE (Freeland et al. 1984) and NCODE (Huggett et al. 1992) experiments. At QE6, the only existing deployment extends from October 1983 through June 1984.

For each model station the kinetic energy has been partitioned into contributions from barotropic and baroclinic modes, according to a procedure outlined in appendix B. The modes are derived from a flat-bottom theory, which is strictly inapplicable over variable topography. They are, however, still useful because they form an orthogonal set of basis functions that provides an efficient decomposition of the model currents in several instances.

Results of the modal partition applied to experiment 14u are provided in Table 5, where contributions from modes 3 and higher have been grouped together. At W02 in northern Hecate Strait, *M*_{2} currents from experiment 14u are similar to those of 14v (shown in Fig. 6) and, as expected, completely dominated by the barotropic component. At other sites considered, the baroclinic contribution is more significant and usually dominated by the first mode.

*M*_{2} currents at two moorings situated in Queen Charlotte Sound, G04 and Q10, are shown in Figs. 13a and 13b. Station G04 is located on the southern flank of Goose Island Bank, in a region subject to a southward radiating flux of baroclinic energy (see Fig 12b). Ellipses from 14v are nearly identical to those presented by BFCJ96 for this site and represent the barotropic response. In the stratified case, the baroclinic and barotropic modes have comparable kinetic energy levels. Model ellipses capture the observed clockwise rotation of phase with depth, a structure expected from superposition of barotropic and first baroclinic mode components (Farmer and Freeland 1983, p. 171). Station Q10, located in the trough between Middle Bank and Goose Island Bank, is exposed to a baroclinic flux radiating shoreward from the shelf break (Fig. 12b). In the stratified case, the vertical structure, dominated by the first baroclinic and barotropic modes, reproduces relatively accurately the observed *M*_{2} ellipses.

Figures 13c and 13d compare current ellipses at two sites in proximity of Cape St. James. Station QB2 is over the continental slope on the west side of the Queen Charlottes, lying in the path of baroclinic fluxes radiating northwestward along the slope from a source region off the Cape (Fig. 12a). Here the response of the stratified model leads to near-surface ellipses that are considerably larger than the unstratified response, consistent with observations. At station QA3, near the shelf break, observed ellipses have a clockwise rotation in orientation and phase with depth, which is captured, albeit with some differences, in experiment 14u.

Results from a site near Dixon Entrance are presented in Fig. 13e. Station QF2 is located near the mouth of the channel, at a point in the path of southward baroclinic fluxes emanating from Learmonth Bank (Fig. 12c). The model ellipses in both cases are dominated by the nearly rectilinear, barotropic component. However, due to a stronger clockwise component, near-surface ellipses from 14u have an increased eccentricity, consistent with observations. Comparisons with observations at station QF1 (not shown), located on the opposite side of Learmonth Bank (Fig. 12c), yield similar results. Generally, the eccentricity of near-surface ellipses throughout much of Dixon Entrance is increased by the presence of the internal tide.

Finally, we show results in Figs. 13f and 13g for two open ocean stations, QE6 located in 1100 m of water off Dixon Entrance and B03, about 80 km southwest of Brooks Peninsula (Fig. 1a), in 2300 m of water. Results from QF4 and B02 (not presented) are qualitatively similar to QE6 and B03, respectively. At QE6 (which lies along the *x*–*z* section described in section 5a), the offshore-directed internal tidal energy (Fig. 12a) produces a large clockwise-rotating ellipse at the depth of the near-surface current meter. By comparison the unstratified response is characterized by relatively weak counterclockwise motion. Similar comments apply to the B03 site, where the internal tide of the model affects near-surface motions, consistent with observations.

Drakopoulos and Marsden (1993) used a ray tracing approach to determine the pattern of characteristics for the B02 and B03 stations. Their results indicated that a beam of internal energy should be detected at a depth of 1000 m for the B03 site, which is inconsistent with both observations and our model. The authors noted that complex local topography and the number of possible generation sites compound the difficulties of applying ray tracing methods in this region.

## 6. Conclusions

This study is concerned with simulation of tides over the northern coastal waters of British Columbia with a three-dimensional primitive equation model. Two basic configurations of the model have been considered, one with homogeneous water and the second with a horizontally homogeneous, vertically stratified fluid, representative of summer conditions. Tidal forcing is introduced along open ocean boundaries with four leading constituents of the barotropic tide derived from a global tidal atlas. Sensitivity to variations in model resolution and to a representative winter stratification have been explored.

The barotropic response of the model is similar to previous tidal simulations of the region, and comparisons with tide gauge data indicate a relatively accurate representation of surface elevation fields. Tidal currents in the unstratified case also agree closely with the earlier studies, but capture little of the vertical structure observed with the semidiurnal constituents in many regions.

With a stratified medium, an internal tide is established which has a significant influence on tidal currents. The rms magnitude of the *M*_{2} baroclinic current near the surface delineates those regions most strongly affected by the internal tide. Most notably, these are Queen Charlotte Sound and Dixon Entrance, where the shelf is broad and exposed to the open ocean. In Hecate Strait, which is sheltered from the open Pacific by the Queen Charlotte Islands, the internal tide has only limited influence. Semidiurnal *M*_{2} tidal ellipses were compared with observations drawn from a series of moorings situated in regions where an internal tide is indicated. The vertical structure of the model ellipses is generally consistent with the observations from a number of widely separated shelf and slope sites. These comparisons suggest that, despite an idealized uniform stratification, the essential characteristics of the internal tide are retained by the model. Further improvements may require a more accurate initialization, taking account of inhomogeneities in the ambient stratification.

Baroclinic energy fluxes computed from the model indicate that semidiurnal baroclinic tides radiate offshore from localized sources along the continental slope. At these generation sites the barotropic tidal ellipses have a pronounced cross-slope component. The integrated offshore flux is about 0.5 GW, which gives an average flux of about 625 W m^{−1} with a model coastline of 800 km. The analytical, two-dimensional model of Baines (1982) provides (to our knowledge) the only estimate of baroclinic energy generation available for comparison. The magnitude of the offshore flux obtained here exceeds this earlier estimate by a factor of at least 5. The discrepancy may be due to difficulties inherent in application of a two-dimensional model to a region with localized sources of baroclinic energy.

Satellite altimetry has permitted a precise determination of the global rate of tidal energy dissipation for the largest constituents (Cartwright and Ray 1991). However, the long standing problem of accounting for these losses in the ocean remains open. In particular, the significance of losses to internal tides is uncertain and disparate estimates have been made. For example, Baines (1982) calculated that internal tide generation along the continental margins accounts for a negligible fraction (0.3%) of the *M*_{2} dissipation rate. Applying the same analytic model to midocean ridges, Morozov (1995) estimated that baroclinic conversion could represent 25% of the global dissipation. The results from this study suggest that incorporation of stratification effects in three-dimensional models may be necessary to simulate in a realistic fashion the dissipation of semidiurnal constituents over continental shelves. For our coastal region we find that stratification leads to significant additional barotropic losses as a consequence of the generation of internal tides.

## Acknowledgments

We thank Mike Foreman for helpful discussions on tidal modelling. We have also benefitted from advice provided by Bill Crawford, Howard Freeland, and Rick Thomson. Contributions to software development were made by Jane Eert, Juijing Gu, and Dave Ramsden. This project was partially supported by the Canadian Panel of Energy Research and Development, Project 67146. We acknowledge support from the Pittsburgh Supercomputing Center.

## REFERENCES

_{2}), atlas of tidal charts and maps. Naval Surface Weapons Center Rep. NSWCTR 79-414.

_{−2}), atlas of tidal charts and maps. Naval Surface Weapons Center Rep. NSWCTR 81-122.

_{−1}), atlas of tidal charts and maps. Naval Surface Weapons Center Rep. NSWCTR 81-142.

_{−1}), atlas of tidal charts and maps. Naval Surface Weapons Center Rep. NSWCTR 81-144.

### APPENDIX A

#### Calculation of Energy Fluxes

The equation governing the rate of change of energy of an incompressible fluid is (Gill 1982; section 4.7)

where *E*_{T} = *ρ*[(1/2)|**u**|^{2} + *gz* + *E*] is the sum of kinetic, potential, and internal energy, with **u** the three-dimensional velocity vector; *G* and *D* represent forcing and dissipative effects, respectively. The energy flux density,

has contributions from the work done by pressure forces, advection of energy, and additional terms involving viscous fluxes and other nonconservative processes. Integration of (A1) over a fixed control volume and time averaging yields a balance between mean fluxes normal to the volume surface and mean forcing and dissipation,

where the angle brackets denote a time-average over period *T,*

We consider the contributions to the mean horizontal energy flux density due to the rate of working by pressure forces associated with barotropic (subscript bt) and baroclinic (subscript bc) tidal motions,

Additional contributions to the flux from advection and diffusion are neglected as these are expected to be of second-order for a small amplitude wave field. For these calculations, the horizontal velocity field **u**_{h} is decomposed into a depth-averaged barotropic component

and an internal baroclinic component

Here **u**_{bc} is interpreted as current associated with the baroclinic tide although it may be contaminated by viscous effects in the bottom boundary layer (see comments below).

The barotropic pressure field *p*_{bt} is associated with displacements of the sea surface, *η*,

A method similar to that of Holloway (1996) has been adopted for calculation of the baroclinic fluxes. Pressure perturbations associated with internal tidal motions are assumed to be governed by hydrostatic balance ∂_{z}*p* = −*ρ**g* and a linearized adiabatic density equation ∂_{t}*ρ* = −*w*∂_{z}*ρ*. The baroclinic pressure field is then obtained from the relation

where *w* is the vertical velocity. A barotropic component, *w*_{bt} = (1 + *σ*) *D**η*/*Dt* + *σ* (**u**_{bt}·**∇***H*) with *σ* = (*z* − *η*)/(*H* + *η*), may be removed from *w,* but this has been verified to make little difference to the baroclinic fluxes.

To evaluate time averages it is useful to expand the velocity and pressure fields in a series of harmonics. For example, consider a rotation through an angle *θ*_{n} to define coordinates (*x*′, *y*′), such that the *x*′ axis is oriented in the direction of the positive semimajor axis of the local tidal ellipse. Then barotropic and baroclinic velocities and the perturbation component of the pressure fields for each constituent are of the form

Here (*u*^{′}_{n}, *υ*^{′}_{n}) are current amplitudes at frequency *ω*_{n} in the (*x*′, *y*′) directions respectively; *g*_{n} is the phase lag of maximum current behind maximum tidal potential; *p*_{n} and *ϕ*_{n} are the pressure amplitude and phase lag respectively; and *θ*_{n} is the orientation of the semimajor axis with respect to (*x*, *y*) coordinates. Using the expansion (A8), it is possible to express components of the rotated mean energy flux vector **M**^{n}_{x′} = 〈**u**′*p*〉_{n} as

A rotation must be applied to recover the energy fluxes in the (*x*, *y*) coordinate of the model:

Terms involving products of velocities and pressures at different frequencies have been omitted because these tend to zero as the averaging period *T* increases. Thus, provided *T* is large compared to the constituent periods and that conditions are stationary, the total mean energy flux is the sum of fluxes at each frequency **M** = Σ_{n}**M**^{n}_{x}. Finally, the depth-integrated energy flux vectors presented in sections 4 and 5 are given by the following expressions, respectively:

The harmonic amplitudes and phases of (A8) and the *θ*_{n} were extracted for each tidal constituent using the least squares fit method of Foreman (1977, 1978). For the barotropic energy flux, these parameters are two-dimensional fields, defined at each horizontal grid point of the model. For the baroclinic fluxes, they are three-dimensional, additionally dependent on the vertical coordinate.

As noted above, two processes contribute to **u**_{bc}: the baroclinic tide, which we wish to isolate, and unwanted viscous effects in the bottom Ekman layer, which should be most significant over shallow shelf regions. To examine whether viscous effects significantly contaminate the energy fluxes, these were recalculated with the baroclinic current in (A3) replaced by (**u**_{h} − **u**_{bt})|_{14u} − (**u**_{h} − **u**_{bt})|_{14v}. That is, the “internal” component of the current from the unstratified experiment (14v), which is due to bottom friction, is subtracted from the baroclinic current of the stratified case (14u). The recalculated depth-integrated fluxes were visually indistinguishable from those given in Figs. 12a–c. The net integrated flux from the shelf into the deep ocean was virtually unchanged. It seems reasonable then to assume that there is no significant contamination of the energy fluxes due to viscous effects.

### APPENDIX B

#### Modal Decomposition

The following procedure was applied to partition the kinetic energy into barotropic and baroclinic modes at the constituent frequencies of the tide. For each station the horizontal currents are expanded as

where *k* is the mode index and *K* the number of vertical grid points. The Π_{k}(*z*) are eigensolutions of the Sturm–Liouville problem,

with boundary conditions Π^{′}_{k} = 0 at *z* = 0, −*H. *Equation (B2) is appropriate to internal wave motions over a flat bottom in the hydrostatic (*ω*^{2} ≪ *N*^{2}), Boussinesq limit (LeBlond and Mysak 1978); *N* is the buoyancy frequency and *λ*^{2}_{k} the eigenvalue.

The modes are normalized according to the condition

where *δ*_{k,j} is the Kronecker delta. Modal currents (*û*_{k}, *υ̂*_{k}) are thus given by

The (*û*_{k}, *v̂*_{k}) are harmonically analyzed in the usual manner and mean kinetic energy in each mode for constituent *n* is obtained as (*A*^{2}_{n,k} + *B*^{2}_{n,k})/4, where *A*_{n,k} and *B*_{n,k} are the semimajor and semiminor axes lengths of the modal current ellipse, respectively.

## Footnotes

*Corresponding author address:* Dr. Patrick F. Cummins, Institute of Ocean Sciences, P.O. Box 6000, 9860 West Saanich Rd. Sidney, BC V8L 4B2, Canada.

Email: pcummins@ios.bc.ca