A Circulation Model for the Discovery Islands, British Columbia

Abstract

A finite volume, ocean circulation model was applied to the Discovery Islands region of British Columbia and used to simulate the three-dimensional velocity, temperature and salinity fields required by a companion biological transport model. The circulation model was initialized with a combination of climatological data and recent temperature and salinity observations, and forced with i) winds measured at seventeen weather stations, ii) the discharges from twelve rivers, and iii) five tidal constituents. A simulation for the period 1 April to 28 April 2010 was evaluated using simultaneous observations from three current meter moorings and the harmonics computed from historical measurements at twenty-four tide gauges. Though the model tidal elevations were shown to be in excellent agreement with the observations, profiles of model tidal speed versus depth generally did not capture observed vertical variations as well. Mean and low-pass filtered flow fields, though reasonably accurate near the surface, were also found to deteriorate farther down the water column. However, near-surface, model current harmonics were shown to be in reasonable agreement with those used to produce annual predictions at five sites in the Canadian Tide and Current Tables. Though the winds were not found to be a significant contributor to mean flow fields over the simulation period, tidal rectification was. Numerous residual and transient eddies that may lead to retentive regions in subsequent Lagrangian studies were predicted by the model. Future work and improvements to overcome model deficiencies are briefly outlined.

RÉSUMÉ [Traduit par la rédaction] Nous avons appliqué un modèle de volume fini des courants océaniques à la région des îles Discovery en Colombie-Britannique et nous nous en sommes servis pour simuler les champs de la vélocité, de la température et de la salinité en trois dimensions, exigés par un modèle de transport biologique d'accompagnement. Le modèle de circulation a été initialisé au moyen d'une combinaison de données climatologiques et d'observations récentes liées à la température et à la salinité et il a été forcé avec les vents mesurés à 17 stations météorologiques, le débit de 12 cours d'eau et cinq composantes de marées. Nous avons évalué une simulation pour la période allant du 1^er avril au 28 avril 2010 à l'aide d'observations effectuées simultanément par trois courantomètres, et les harmoniques ont été calculées à partir des mesures effectuées dans le passé à 24 marégraphes. Bien que les hauteurs de la marée établies par le modèle concordent parfaitement avec les observations, les profils de la vitesse de la marée dans le modèle en fonction de la profondeur ne réflétaient généralement pas aussi bien les variations verticales observées. De plus, nous avons constaté que les champs de courant moyens et auxquels on avait appliqué un filtre passe-bas perdaient en précision plus bas dans la colonne d'eau bien que leur valeur près de la surface soit raisonnablement exacte. Cependant, les harmoniques des courants près de la surface, prédites par le modèle concordaient relativement bien avec celles qui avaient été utilisées dans le but d’établir des prévisions annuelles aux cinq emplacements figurant dans les Tables des marées et courants du Canada. Nous avons constaté que les vents ne représentaient pas un apport important par rapport aux champs de courant moyens au cours de la période de simulation, à la différence de la rectification maréale. Le modèle a prédit de nombreuses perturbations transitoires et résiduelles, qui pourraient permettre d’établir les régions de rétention dans les études ultérieures de type Lagrangien. Nous décrivons brièvement les travaux et les améliorations à faire pour corriger les lacunes du modèle.

Keywords:

• circulation model
• Discovery Islands
• FVCOM
• April 2010 simulation

1 Introduction

The Discovery Islands region, between east-central Vancouver Island and the British Columbia mainland (Fig. 1) of Canada, is a complex network of narrow channels and deep fjords whose currents comprise wind-driven flows, estuarine flows largely driven by seasonal freshwater discharge, and some of the strongest tidal currents in the world (Canadian Tide and Current Tables, 2010). The Discovery Passage and Johnstone Strait corridor is a major tugboat, barge and cruise ship route between Alaska and ports such as Vancouver and Seattle in the Salish Sea, and many waterways within the Discovery Islands region are favourite destinations for pleasure craft seeking to enjoy the scenery and generally abundant fishing. Numerous salmon farms are also found in the region and, as in the Broughton Archipelago to the northwest, their presence has given rise to environmental concerns (Krkošek et al., 2005, 2007, 2008; Brooks and Stucchi, 2006; Brooks and Jones, 2008) and a need to understand better the circulation patterns controlling the dispersion of organisms and pollutants. Two particular farm-related issues that are of immediate concern are the dispersion of the infectious haematopoietic necrosis virus (IHNV) between farms (Saksida, 2006) and the possibility that sea lice originating on the farms may be adversely affecting juvenile salmon migrating from the Fraser (Peterman et al., 2010) and other rivers (Price et al., 2010) to the Pacific Ocean.

Fig. 1 Map of the Discovery Islands region showing relevant place names, rivers, bodies of water, weather stations, CTD sites, Tide Table sites and the Barnes Bay fish farm.

The physical oceanography of the Discovery Islands region has not been widely studied. Though tide gauges have been maintained at numerous locations for varying lengths of time (e.g., Fig. 2 and Table 1), moored current meter measurements and salinity and temperature profiles have been sparse in the central portion of the region, largely because of the strong tidal currents and the potential danger their deployment and recovery pose to small vessels. Based on tide, current and water property observations taken in Johnstone Strait, at the northern extent of the Discovery Islands region, Thomson and Huggett (1980) found westward mean flows with magnitudes up to 30 cm s⁻¹ in the upper 100 m of the channel; eastward mean flows with magnitudes up to 20 cm s⁻¹ in the lower 200 m; M₂ along-channel baroclinic current amplitudes as large as 50 cm s⁻¹; and a relatively high (>4.5 mL L⁻¹) and near-uniform distribution of dissolved oxygen observed throughout the year, suggesting vigorous tidal pumping and turbulent mixing along the entire length of the channel. In 1987–88, Woodward and colleagues combined short-term (typically two to three weeks) shipboard current observations with longer bottom pressure measurements in many narrow passes in order to estimate the tidal harmonics that would be the basis of predictions for the Canadian Tide and Current Tables (2010). However, there were problems in assigning these harmonics to specific locations, and this will be discussed later when the model is evaluated.

Table 1. a) M₂ and b) K₁ model tidal elevation amplitudes (cm) and phases (degrees, gmt) versus values arising from tide gauge analyses. Length (days) is the tide gauge time series length while D (cm) is the difference between modelled and observed values in complex space. Tide gauge numbers are consistent with Fig 2.

No.	site	length	Tide gauge		Model		D
			amplitude	phase	amplitude	phase
a)
1	Little River	365	99.6	33.9	99.2	33.7	0.6
2	Powell River	369	100.7	34.3	100.4	33.7	1.0
3	Lund	224	102.2	35.4	100.8	35.1	1.4
4	Twin Islands	369	101.3	35.4	99.5	35.0	1.9
5	Campbell River	365	76.2	14.4	75.6	20.2	7.7
6	Nymphe Cove	539	61.5	350.4	57.3	351.8	4.4
7	Brown Bay	365	93.5	315.9	91.9	315.5	1.7
8	Seymour Narrows	369	94.6	320.1	95.7	319.6	1.4
9	Maude Island East	368	55.6	7.4	52.1	13.9	7.0
10	Welsford Island	364	99.4	35.1	98.2	38.1	5.4
11	Redonda Bay	339	97.5	36.7	99.4	36.3	2.0
12	Channel Islands	288	102.6	35.9	101.2	36.2	1.5
13	Chatham Point	369	90.3	305.1	90.0	305.0	0.3
14	Okis Islands	366	92.0	312.6	89.6	313.0	2.5
15	Owen Bay	369	85.0	319.9	83.4	319.8	1.6
16	Bodega Anchorage	365	72.9	5.7	72.0	348.7	21.4
17	Big Bay	369	75.5	14.9	71.9	10.6	6.5
18	Turnback Point	365	102.0	37.0	99.4	37.5	2.8
19	Orford Bay	282	101.5	37.2	99.9	37.2	1.5
20	Waddington Harbour	169	103.4	38.0	100.8	37.5	2.7
21	Shoal Bay	387	89.8	307.5	88.2	307.3	1.6
22	Cordero Islands	386	85.8	304.4	85.6	302.4	3.0
23	Kelsey Bay	375	117.0	276.3	112.6	274.3	5.9
24	Yorke Island	362	117.1	271.8	112.1	271.2	5.1
b)
1	Little River	365	90.5	287.5	87.6	288.4	3.3
2	Powell River	369	90.4	286.6	87.3	287.3	3.3
3	Lund	224	88.9	287.9	88.2	288.1	0.8
4	Twin Islands	369	90.4	287.5	87.9	288.3	2.8
5	Campbell River	365	82.6	283.7	82.8	285.8	3.1
6	Nymphe Cove	539	77.0	279.9	75.3	282.4	3.8
7	Brown Bay	365	67.9	270.1	66.5	270.4	1.5
8	Seymour Narrows	369	69.2	272.1	69.1	271.2	1.2
9	Maude Island East	368	81.1	283.9	78.9	286.2	3.9
10	Welsford Island	364	91.1	286.9	88.9	289.1	4.1
11	Redonda Bay	339	87.1	287.4	88.7	288.6	2.4
12	Channel Islands	288	89.9	288.0	88.4	288.6	1.7
13	Chatham Point	369	65.4	270.5	64.5	270.2	1.0
14	Okis Islands	366	67.5	269.8	66.1	271.1	2.0
15	Owen Bay	369	67.8	272.7	67.4	273.3	0.8
16	Bodega Anchorage	365	79.8	281.9	73.1	280.5	7.0
17	Big Bay	369	83.3	283.5	79.5	283.8	3.8
18	Turnback Point	365	91.7	287.6	88.6	288.9	3.7
19	Orford Bay	282	90.3	288.1	88.8	289.0	2.1
20	Waddington Harbour	169	89.2	288.2	88.9	288.9	1.2
21	Shoal Bay	387	66.6	269.6	65.1	270.7	2.0
22	Cordero Islands	386	66.1	269.9	64.2	270.7	2.2
23	Kelsey Bay	375	57.7	261.4	56.7	261.6	1.0
24	Yorke Island	362	55.7	260.0	55.7	261.9	1.9

Fig. 2 Locations of tide gauges (numbered circles) and current meters (red triangles) used to evaluate model accuracy. Colour contours are bathymetry (m). Current meters are Nodales Channel (NC), Discovery Passage (DP) and Cape Mudge (CM).

Numerical modelling in the Discovery Islands region has also been sparse. Though Crean et al. (1988) employed a network of single grid cells to represent the three major channels exiting the northern end of the Strait of Georgia (SOG), the major focus of their studies was the Georgia Strait, Juan de Fuca Strait, Puget Sound system, and their model was not tuned to capture the Discovery Islands flows accurately. Foreman et al. (2004) included the Discovery Islands region in their data assimilating, barotropic, finite element model of the waters around Vancouver Island, but in this case the primary focus was M₂ tidal dissipation. Sutherland et al. (2007) used a similar two-dimensional (2-D) model to estimate the tidal power potential obtained by inserting turbine farms in Johnstone Strait, Discovery Passage and Cordero Channel and found them to be very close to the analytical estimates provided by Garrett and Cummins (2005). This was a refinement of similar studies carried out for BC Hydro (2002) for all British Columbia waters. Finally Jiang and Fissel (2010) used a much more finely resolved grid and a full three-diemnsional (3-D), baroclinic finite difference model to assess potential turbine locations in southern Discovery Passage, including the possibility of partially or fully removing a dam in Canoe Pass.

The modelling carried out in this study employs the Finite Volume Coastal Ocean Model (FVCOM) developed by Chen et al. (2003, 2004, 2006a, 2006b). FVCOM solves the 3-D hydrodynamic equations for velocity and surface elevation and 3-D transport-diffusion equations for salinity and temperature in the presence of turbulent mixing. It also employs a variable resolution triangular grid to provide much more flexibility than a regular rectangular grid in representing complicated regions such as the coastline and bathymetry in the Discovery Islands region. It can be forced with any combination of tides, wind, river runoff, surface heating and open ocean inflows, and though it has coupled nutrient-phytoplankton-zooplankton and sediment transport components, they are not employed here. This is the same code that Foreman et al. (2009) used to simulate the circulation in the Broughton Archipelago region, just to the north of the Discovery Islands and that was coupled to a sea lice behaviour and dispersion model (Stucchi et al., 2011). FVCOM was also applied to simulate 3-D baroclinic circulation and vertical stratification in Placentia Bay off the east coast of Canada (Ma et al., 2012). Further details on FVCOM and a more complete list of publications that have used it can be found at MEDM (2004).

The time period for the FVCOM simulation described herein was 1–28 April 2010. This period coincides with physical sampling programs that allowed i) the estimation of reasonably accurate 3-D salinity and temperature initial conditions, ii) the use of direct wind observations as model surface forcing, and iii) the evaluation of accuracy using acoustic Doppler current profiler (ADCP) measurements from three moorings in the region. Running under Message Passing Interface (MPI), each simulation for this 28-day period required approximately 38 hours on 64 processing units of an Intel cluster machine.

This paper is organized as follows. In Section 2, the model setup; tidal, river discharge and wind forcing; and the calculation of initial salinity and temperature fields are briefly described. In Section 3, the results of the FVCOM simulation are presented and evaluated using both simultaneous and previous observations. Finally, in the last section, the results are summarized and future work is outlined.

2 Model setup and forcing

a Model Geometry

The triangular grid (Fig. 3) for this FVCOM application has 35609 nodes, 65473 elements and a horizontal grid resolution that varies from approximately 90 m in Seymour Narrows to 1.7 km in the SOG. In order to allow the freshwater discharges of twelve rivers to adjust dynamically before they reach the ocean, generic channels of approximately 5 km length were included in the grid. (Further details on prescribing these discharges will be given in Section 2c.) In the vertical, a sigma-coordinate system was used with 21 levels satisfying a hyperbolic tangent distribution (Pietrzak et al., 2002) with high resolution at the surface and a maximum interlayer spacing of 15% of the depth at the bottom. The grid was generated with an update of the TriGrid software developed by Henry and Walters (1993). The model bathymetry was generated from Canadian Hydrographic Service single-beam digital charts with some sections updated with values from multi-beam surveys. The depths ranged from 5 m in coastal regions and rivers to over 700 m within Bute Inlet. While the true depth in some coastal sections approached zero, specifying a minimum depth of 5 m economized on the model run times by avoiding the creation of “mud flats” in shallow areas near islands and estuaries and allowing freshwater emanating from the heads of rivers to mix dynamically before reaching the river mouths. Though FVCOM has a provision for wetting-and-drying, the generally steep-sided nature of shorelines in the region (parts of Bute Inlet are over 700 m deep and carved into mountainous terrain rising up to 3000 m) means that, unlike the low-lying estuarine application described in Zhao et al. (2010), mudflats play a very minor role in the overall circulation dynamics and can be ignored.

Fig. 3 Horizontal grid and rivers whose discharges were used for the FVCOM simulations.

The model bathymetry was smoothed with a volume preserving technique (Mellor et al., 1994) that limits the ratio Δh/h within each triangle, where Δh is the difference between maximum and minimum depth and h is the average depth. The intent is to avoid “hydrostatic inconsistency” (Haney, 1991), a potential problem with sigma-coordinate (or terrain-following) oceanographic models in the presence of steep bathymetry. If σ denotes the fractional vertical spacing (−1 ≤ σ ≤ 0), Sikirić et al. (2009) recommended using the regional oceanic modelling system (ROMS; Shchepetkin and McWilliams, 2005) grids with a “safe” hydrostatic inconsistency (or Haney) number, rx ₁= (Δh/h)(σ/Δσ), less than 3 but stated that most simulations are carried out with rx ₁ ≤ 6 and Δh/h ≤ 0.4. Though ROMS findings are not necessarily transferable to FVCOM applications, trial-and-error simulations and discussions with Tarang Khangoandar (personal communication, 2011) indicated that a sigma-coordinate choice with (σ/Δσ) ≤ 0.16 and (Δh/h) ≤ 0.6 would retain relatively accurate bathymetry and inhibit excessive mixing. Consequently, those were the criteria used here. Apart from the river channels, this smoothing resulted in most coastal depths becoming larger than the 5 m minimum.

Figure 4a shows original and smoothed bathymetry along a Discovery Passage cross-channel transect that includes mooring DP in Discovery Channel (Fig. 2) while Fig. 4b shows a similar transect extending from near the mouth of the Homathko River channel into the upper reaches of Bute Inlet. In both cases the Δh/h < 0.6 smoothing is seen to capture the original bathymetry more accurately than Δh/h < 0.3 smoothing and is not much worse than with Δh/h < 0.9. A similar smoothing value (0.4) was used for the Placentia Bay model (Ma et al., 2012).

Fig. 4 Original and smoothed bathymetry (m) for a) a Discovery Channel transect through mooring DP and b) an along-channel transect from the lower reaches of the Homathko River into Bute Inlet. Δh/h values denote smoothing criteria.

b Numerical Considerations

Numerical integrations in FVCOM are carried out with a second-order accurate finite volume flux discrete scheme where internal and external modes are “split” and integrated over distinct time steps. The version of FVCOM employed in our study (2.7.1) had explicit time stepping, and trial-and-error testing determined that an external time step of 0.075 s preserved stability over our simulation. The time step for the internal mode is calculated as Δt_I  = Δt_EI _SPLIT with I _SPLIT = 10 in our application.

A Smagorinsky eddy parameterization with coefficient C = 0.2 was used to provide horizontal turbulence closure. The variation of q-ϵ turbulence closure that Tian and Chen (2006) found to be most accurate for Georges Bank and that Foreman et al. (unpublished manuscript) also found to be best for the highly stratified waters of the Broughton Archipelago was also employed here. The background vertical diffusion and viscosity was set to 10⁻⁶ m² s⁻¹, and the model bottom drag coefficient was expressed using a logarithmic wall-layer formulation with the bottom roughness parameter set to 0.001 and a minimal value of 0.0025.

c Initial Conditions and Lateral Boundary Forcing

Initial three-dimensional fields of salinity and temperature were computed by combining seasonal historical climatology (Foreman et al., 2008) interpolated to 1 April with conductivity, temperature, depth (CTD) measurements taken over the period 1–4 April 2010. The locations of these CTD stations are shown in Fig. 1. The climatological values were first interpolated to the 3-D model grid using a Kriging method within the Surfer modelling package (Golden Software, Inc., 2006), then the more recent values were incorporated into these fields with an irregular-grid diffusion solver. Figure 5 shows the resultant surface salinity (practical salinity scale used) and temperature fields. The temperature variation is seen to be about 2.5°C with generally colder values at the heads of the fjords. The salinities also tend to be lower (fresher) in the fjords, largely because of river runoff, and higher in regions with strong tidal mixing.

Fig. 5 Initial surface salinity and temperature (°C) interpolated from historical climatology and recent observations at the CTD locations shown in Fig 1.

Although only the results from tidal constituents M₂ and K₁ are reported here, the five constituents M₂, S₂, N₂, K₁ and O₁ were included in all model simulations. This was done to ensure that the major non-linear interactions, primarily via quadratic bottom friction and advection, would be approximately correct. Though overtides like M₄, MS₄ and M₆ were generated by FVCOM and found to be relatively large in regions like Discovery Passage, they will not be described here. At Campbell River (station #5 in Fig. 2), M₂, S₂, N₂, K₁ and O₁ account for 74% of the tidal height range in the diurnal and semi-diurnal frequency bands. Boundary conditions for the tidal and mean components of the circulation model are zero flow normal to the coast and specified elevations along the Johnstone Strait, SOG and Malaspina Strait open boundaries. Values along these boundaries were approximated using a combination of tide gauge harmonics used to produce the Canadian Tide and Current Tables (2010) and results from larger domain models (Foreman et al., 2004).

Though baroclinic waves and internal pressure gradients can be expected to evolve temporally along the open boundaries, the absence of both specific observations for our simulation period and output from a model whose domain encompasses this one, necessitated the implementation of a condition that nudged the salinities and temperatures along these boundaries back to their initial conditions. Though admittedly not correct, because these boundaries are relatively far from the central Discovery Islands region of interest and this region has strong tidal mixing, the impact of missing, temporally evolving, baroclinic effects along the boundaries should be minimal over our 28-day simulation period.

d Atmospheric Forcing

As with the Broughton Archipelago model domain (Foreman et al., 2009), there are few permanent weather stations in the Discovery Islands region and no atmospheric models with sufficient resolution to capture the expected wind steering accurately along narrow passages. Twelve weather stations were deployed between October 2009 and February 2010 to augment the five already in place (Fig. 1), and they have been maintained with few problems. Observed hourly winds for 1–28 April at three stations close to the current meter moorings (Fig. 6) show considerable variation over a relatively small region. The Sentry Shoal winds in the northern SOG are seen to be generally northwestward up to 8 April, southeastward for the next nine days, and then variable for the remainder of the month. The Cinque Island winds near current meter mooring DP (Fig. 2) generally display the same patterns, variability and strength but have directions that show stronger alignment with Discovery Channel. The fact that the Cinque Island wind events appear more persistent simply reflects that the Sentry Shoal time series has more missing data. However, the Lee Islands winds near the current meter mooring in Nodales Channel (NC) are seen to be much more variable and generally weaker; they also display little similarity to the wind patterns at the other two locations. This time series also has more missing data than the Cinque Islands time series. The measured winds from all the weather stations were interpolated or extrapolated to model grid elements with a thin plate spline technique (Wahba, 1990) available from MATLAB (2009). Because many of the channels in the Discovery Islands region are not as steep-sided as, and thus do not steer the winds to the same degree as, those in the Broughton Archipelago, this approach, rather than the representer interpolation technique employed in Foreman et al. (2009), was used here.

Fig. 6 Winds on 1–28 April measured at the Sentry Shoal, Cinque Islands and Lee Islands weather stations shown in Fig 1.

Only one long time series of temperatures is available for the Discovery Islands region. Mainstream Canada measured water temperatures once a day at depths of 1, 5, 10 and 15 m at their Barnes Bay aquaculture site (Fig. 1) from 1 June 2009 to 30 June 2010 (Fig. 7). For the 1–28 April period of our model simulation, the conditions were mainly rainy or overcast; there was little solar radiation, and the air temperature remained relatively close to the sea surface temperature which ranged between 8.0° and 8.7°C. Heat flux forcing was, therefore, assumed to be unimportant and was neglected in the model simulations. However, it should be included in subsequent summer simulations because Fig. 7 shows the Barnes Bay temperatures rising to 12°C and having oscillations of 1.5°C over 10-day periods in August 2009.

Fig. 7 Daily water temperatures taken at 1 and 15 m depth at the Barnes Bay aquaculture site for the period 1 June 2009 to 30 June 2010. The grey strip denotes the 1–28 April model simulation period.

e Freshwater Forcing

The relatively strong estuarine circulation that Thomson and Huggett (1980) found in Johnstone Strait suggests significant sources of fresh water upstream, in either the Discovery Islands region or the SOG or both. Though strong tidal mixing in most of the channels between the northern SOG and Johnstone Strait can be expected to diminish the stratification and estuarine flows that large river discharges such as those generated by the Fraser River in the SOG, the extent to which the Johnstone Strait estuarine flows arise from stratification in the SOG or river discharges in the Discovery Islands region is not clear. To model these estuarine systems accurately, daily freshwater discharge was specified for nine rivers on the British Columbia mainland side of the model domain: the Powell, Toba, Brem, Southgate, Homathko, Estero, Phillips, Apple and Stafford rivers; along with the Salmon, Campbell and Oyster rivers on the eastern coast of Vancouver Island (Fig. 3).

All the freshwater sources were incorporated into the FVCOM model as river discharges divided equally across four nodes separated by approximately 50 m at the heads of each of the respective river channels and over the minimum depth of 5 m in the vertical. As described in Foreman et al. (2009), this strategy allowed some dynamical adjustment in the channels and produced better results than treating the fresh water as point sources entering abruptly at the coastal boundary. The river discharge data were obtained from the Hydrometric Network (EC, 2010a) maintained by the Water Survey of Canada (WSC; a branch of Environment Canada) for the Homathko, Salmon and Oyster rivers while discharges for the Toba, Brem, Southgate, Estero, Phillips, Apple and Stafford rivers were estimated by comparing watershed areas to that of the Homathko (Table 2). Flows in the Campbell River are regulated by several dams; the John Hart is the closest to the river mouth. Flows from this dam are typically maintained at 124 m³ s⁻¹ (BCRP, 2000) while flows for the Quinsam River, which enters the Campbell River below this dam, are gauged. So the total discharge entering Discovery Passage is estimated as the sum of these two components. Figure 8 shows the Homathko, Salmon, Campbell and Oyster river discharges for 1–28 April 2010. While the Salmon River discharge is highly variable because of rain events, the Homathko River shows more steadily increasing values toward month-end as air temperatures increase and the snowpack begins to melt. The prescribed river temperatures associated with these discharges were estimated from surface CTD measurements taken near the river mouths while, with the exception of the Estero Basin (see discussions in the next section), the prescribed river salinities were all set to zero.

Fig. 8 Discharges for the Homathko, Salmon, Oyster and Campbell rivers (EC, 2010b).

Table 2. Estimated and gauged river discharges for 1 April 2010; e = estimated, g = gauged.

River	Estimated or gauged	Watershed area as fraction of Homathko	Discharge (m^3 s^−1)
Powell	e	0.10	8.0
Toba	e	0.310	24.8
Brem	e	0.046	3.7
Southgate	e	0.385	30.8
Homathko	g	1.0	80.0
Estero	e	0.025	2.0
Phillips	e	0.05	4.0
Apple	e	0.035	2.8
Stafford	e	0.045	3.6
Salmon	g	—	80.0
Campbell	g	—	140.0
Oyster	g	—	20.0

3 Model results

The model was run for 28 days with the wind and boundary tidal forcing ramped up from zero over the first two days. Elevation, velocity, temperature and salinity values were stored at all nodes or elements and layers every hour throughout the simulation, and a harmonic analysis was performed over the last 25 days to extract average fields and tidal energy at constituent frequencies. Comparisons of model results with elevation and current amplitudes and phases computed at the tide gauge and current meter locations shown in Fig. 2 are summarized in Table 1. Elevation differences are calculated as distances, D, in the complex plane; that is,

where A_o , A_m , g_o and g_m are the observed and modelled amplitudes and phases of a particular constituent, respectively. Models with D values within 5% of A_o are generally considered to have acceptable accuracy.

Table 1 gives M₂ and K₁ modelled and observed elevation harmonics at the tide gauge locations shown in Fig. 2. Average D values over all 24 sites are 3.8, 2.5, 1.4, 1.5 and 1.3 cm for M_2, K_1, S_2, O₁ and N₂, respectively. Co-amplitude plots for M₂ and K₁ are shown in Fig. 9. The largest M₂ discrepancies are generally found in regions with large amplitude gradients, such as near Bodega Anchorage, and indicate that the model minima are slightly displaced. Notice that the M₂ phases (Table 1a) in Johnstone Strait (gauge numbers 23–24) are approximately eight hours (240°) larger than those along the southern SOG boundary (gauge numbers 1–4). It is these out-of-phase differences for all the semi-diurnal constituents that cause the large currents in many of the channels in the Discovery Islands region. In essence, these phase differences arise from the different travel times when one part of the M₂ wave propagating northward along the Oregon–Washington coast enters Juan de Fuca Strait and continues up the SOG and meets a second part that has continued along the west coast of Vancouver Island, turned around the northern tip, and propagated into Queen Charlotte Strait and Johnstone Strait. The destructive interference between these out-of-phase northward and southward waves leads to the minimal amplitudes seen in Fig. 9a. The rapid change in M₂ amplitudes adjacent to the Johnstone Strait boundary reflects the fact that the forcing values on this boundary had to be artificially adjusted to produce more accurate values farther to the south and east. The need for this adjustment may, in part, be caused by the sponge layer that was required to maintain numerical stability in the region. A similar adjustment was not needed with K₁, or along the southern forcing boundary.

Fig. 9 Model co-amplitudes (cm) for a) M₂ and b) K₁.

Figure 10 shows model mean flows and surface elevations at 1 m depth computed by harmonic analyses for the period 4–28 April 2010. The patterns are a combination of surface estuarine flows largely driven by the freshwater discharges emanating from the rivers (Fig. 3) and tidally rectified flows driven by strong gradients in the velocity fields. Because a similar plot for a model simulation without wind forcing (not shown) is essentially identical, we conclude that, at least for this time period, the winds have made little contribution to these mean velocities. Only vectors at nodes separated by a minimum of 600 m are shown in the larger region, and though the numerous crossing vectors suggest possible inconsistencies in the flow field, zoom-ins such as that displayed for the Seymour Narrows region (Fig. 1) show that these inconsistencies arise from a rich eddy field. The extent to which these Eulerian eddies retain passive particles, such as viruses emanating from salmon farms and processing plants in the region, will be investigated and reported in future studies. The flow field at 100 m depth (not shown) has a somewhat weakened eddy field and bottom estuarine flows (e.g., in Nodales Channel) in the opposite direction to those in Fig. 10.

Fig. 10 Mean model surface elevations and flows at 1 m depth computed by harmonic analyses for the period 4–28 April 2010. Only vectors at nodes separated by a minimum of 600 m are shown in the larger region; all vectors are shown in the Seymour Narrows insert.

Figure 11 shows the maximum surface speeds over the 28-day model simulation period. Large values in regions such as Seymour Narrows and Arran Rapids (Fig. 1) are primarily caused by the tides, but estuarine flows and river discharges do contribute in regions such as Nodales Channel, Johnstone Strait, and near the mouths of rivers like the Oyster.

Fig. 11 Maximum surface speeds (m s⁻¹) over the 28-day model simulation period.

Figure 12 compares modelled and observed vertical profiles of M₂ and K₁ tidal current amplitude at the three mooring locations shown in Fig. 2. As in all cases the tidal ellipses have very small semi-minor axes, and there is close agreement between the modelled and observed angles of inclination; the flows are essentially rectilinear, and our comparison is restricted to amplitude (maximum speed). Generally, the modelled K₁ amplitudes are seen to be reasonably accurate near the surface and to degrade with depth. While the M₂ model surface speeds are also reasonably accurate at DP (Fig. 12b) they are too large at NC (Fig. 12a), too small at the mooring at Cape Mudge (CM; Fig. 12c), and in general none of the three model profiles capture the observed baroclinicity very well. A run with the Mellor–Yamada level 2.5 (MY2.5) turbulence closure (Mellor and Yamada, 1982) was made to determine the extent to which the mixing parameterization might be influencing these profiles, but it showed little difference from those obtained using the Tian and Chen (2006) scheme. Other possible causes and potential cures for these profile inaccuracies will be discussed in the final section.

Fig. 12 Observed and modelled vertical profiles of M₂ and K₁ tidal speed (semi-major axis, cm s⁻¹) at the current meter moorings (Fig. 2) in a) Nodales Channel (NC), b) Discovery Channel (DP) and c) Cape Mudge (CM).

Model flows in Nodales Channel were found to be highly sensitive to the freshwater discharges prescribed in Estero Basin (Fig. 3). Initial runs with an Estero discharge time series consistent with the values given in Table 2 produced near-surface estuarine and M₂ tidal flows that were much too large at mooring NC (Fig. 2). Because reducing the discharges to almost zero improved matters only slightly, two further sensitivity runs were carried out. In the first, the Estero Basin was removed completely as a source of fresh water while in the second it was retained with half the original discharges, and the incoming salinity was reset to 24.7, rather than the zero value used initially. (The value of 24.7 is the average of the surface values shown in Fig. 5a and those at 5 m depth.) Because these two runs produced similar results that had better near-surface agreement with the NC observations and the second run had more realistic forcing, values from the second run are shown in Figs 10 –13.

Fig. 13 Observed and modelled mean, “along-channel”, vertical profiles at current meter moorings a) NC, b) DP and c) CM. For CM, “along-channel” was chosen as 160° counterclockwise from east, the angle of inclination of the M₂ major semi-axis.

Collectively, these results suggest that it was not so much the Estero Basin discharge that caused the initially inaccurate model flows at NC, rather it was the prescribed salinity. Salinity in FVCOM is computed using a typical advection-diffusion equation. As seen in Fig. 10, the flows near Estero Basin and in Nodales Channel are relatively small compared to elsewhere in the model domain. Consequently, diffusion plays a proportionately larger role in changing salinities from their initial values. Our sensitivity runs suggest that significant estuarine flows can result with no river discharge as long as the prescribed river salinity is substantially different from the ambient values of the inlet into which the river is emptying. Though a similar effect can be expected from prescribed river temperatures, it was avoided in this case by choosing values that were the same as those measured by CTDs near the river mouths and used in computing the initial conditions. A disadvantage of the second sensitivity run is that it introduces a source term to the overall model salt budget, thus violating salt conservation. However, the nudging conditions imposed at the Georgia Strait, Malaspina Strait and Johnstone Strait open boundaries also contribute to the salt (and temperature) budget(s); therefore, precise conservation is not maintained even for the first sensitivity run.

Figure 13 displays modelled and observed vertical profiles for the mean current for the period 4–28 April at the three mooring locations. Though the surface flows at NC (Fig. 13a) are too large and those in the remainder of the southwestward, top-layer, estuarine flow are too small, values below 100 m are reasonably accurate. The thickness of the surface model estuarine layer is seen to be about 70 m versus 90 m for the observed thickness at DP (Fig. 13b); the model flows are too strong virtually all the way down the water column and have a much too weak and thin bottom return flow. Finally, the model mean flow profiles at CM (Fig. 13c) show very little agreement with the observations. This mooring is south of Discovery Passage and its flow directions are not restricted by a relatively narrow channel. Whereas the observed mean flow directions range from northward at 26 m depth to southeastward at 122 m, the model mean flows are westward to northwestward over the entire water column. This disagreement is puzzling because the angles of inclination for the modelled and observed M₂ current ellipses agree very well, both being consistently in the range of 154° to 170° counterclockwise from east. The surface residual circulation computed by Crean et al. (1988, their Fig. 13.15) also showed westward flows around this mooring. Clearly, more investigation is needed to unravel the causes of these disagreements, and this will be discussed in the final section.

Lest one conclude that the sub-tidal flows at these three moorings are relatively constant, Fig. 14 shows observed and modelled low-pass filtered, along-channel currents as functions of depth and time. The ADCP measurements at DP (Fig. 14a, top panels) show a clear spring-neap modulation of the bottom estuarine currents in both magnitude and depth. The associated modelled currents have a similar modulation, though there are clear disagreements in the magnitude, depth range and timing. Spring-neap modulations are not as evident in the NC low-pass filtered currents (Fig. 14b). In this case the model completely misses the estuarine reversal on 7–9 April, perhaps because of gaps in wind observations from the Lee Islands (Fig. 6) weather station. And though Fig. 13c indicates strong disagreement between the modelled and observed mean flows at mooring CM, Fig. 14c shows that there are times when the agreement is reasonable. However, the low-pass filtered ADCP values do exhibit considerable temporal variations while the associated model values are much more persistent.

Fig. 14 Observed and modelled low-pass filtered, “along-channel,” currents (m s⁻¹) versus depth for a) DP (positive is 116° counterclockwise from east), b) NC (positive is 57° counterclockwise from east) and c) CM (positive is 155° counterclockwise from east).

Disagreements between observed and modelled current profiles may be a result of the model inaccurately capturing the vertical stratification. Unfortunately, this is difficult to confirm because, except for the CTD values used to establish the initial temperature and salinity conditions, there were no other observations that could be used to evaluate stratification changes over the period of the simulation. Nevertheless, cross-channel salinity profile changes along transects passing through moorings DP and NC (Fig. 15) suggest that the model stratification has been degraded over the 28-day simulation. Both the near-surface freshwater lenses and the bottom saltier water have been diffused. As mentioned previously, in addition to results from the Tian and Chen (2006) mixing scheme shown here, a simulation was also carried out with MY2.5 (Mellor and Yamada, 1982) turbulent mixing. The fact that its stratification changes showed very little difference from those shown in Fig. 15 suggests that our current shear inaccuracies are not a consequence of the particular general ocean turbulence model (GOTM; Burchard et al., 1999) scheme that was chosen.

Fig. 15 Model salinities for transects crossing Discovery Channel through mooring DP at a) hour 0 and b) hour 672; and Nodales Channel through mooring NC at c) hour 0 and d) hour 672.

Finally, Table 3 compares modelled and observed along-channel surface velocity amplitudes and phases for five locations included in the Canadian Tide and Current Tables (2010; see Fig. 1 for locations). As mentioned earlier, the observed values were not taken with conventional current meters but rather were obtained by combining bottom pressure measurements with shipboard measurements taken from a Zodiac that attempted to maintain a steady position within the flow fields for periods up to eight hours. Though in some cases the Zodiac moved locations depending on whether the flows were in the ebb or flood portion of the tidal cycle, single sets of harmonic constants were produced for each Tide Table location. The precise location and accuracy of these observed harmonics are difficult to estimate; nonetheless, it is interesting to compare them with their model counterparts. Not only does Table 3 show that there is generally good agreement, but it was found that in all cases a slight shift in the observation location would improve that agreement. Except at Hole-in-the-Wall, all model amplitudes are smaller than the Tide Table values while most of the phases are quite close.

Table 3. M₂ and K₁ along-channel model velocity amplitudes (cm s⁻¹) and phases (degrees, gmt) versus values used in the Canadian Tide and Current Tables (2010). Locations are shown in Fig 1.

	M2				K1
	Observations		Model		Observations		Model
	amplitude	phase	amplitude	phase	amplitude	phase	amplitude	phase
Arran Rapids	456	267	409	266	99	149	78	152
Beazley Passage	369	85	305	89	78	321	62	327
Gillard Passage	356	91	270	93	82	333	60	334
Hole-in-the-Wall	373	269	407	269	79	142	82	148
Seymour Narrows	466	113	425	111	101	337	99	345

4 Summary, discussion and future work

We have described the application of FVCOM to the Discovery Islands region of British Columbia and the evaluation of a simulation over the period 1–28 April 2010 using current and sea surface elevation observations. The primary objective for this model development was the production of credible 3-D fields that could be used to study the dispersion of IHNV from aquaculture sites within this region (Saksida, 2006). The model evaluations revealed generally high accuracy in reproducing tidal elevations throughout the region, and while many of the 3-D observed currents were not captured accurately, the modelled near-surface values were generally acceptable. In an attempt to lay the groundwork for future model improvements and simulations, the following discussion highlights the strengths and weaknesses of this FVCOM application and provides some links to the coupled viral transport model that will be described in future manuscripts.

As illustrated in Figs 12, 13 and 14 the primary deficiency in the model simulations was the poor representation of model currents down the water column. These inaccuracies could arise for a variety of reasons, foremost of which are probably the mixing and resolution of the grid and bathymetry. Clearly, this is a very complicated region in terms of coastal irregularities, bathymetry and strong gradients in the flow regimes (Fig. 10). A higher resolution grid that captures the details of these features better needs to be constructed and these simulations repeated. Higher spatial resolution, particularly in regions with large depth gradients, will also reduce inaccuracies arising from the bathymetric smoothing that was needed to avoid hydrostatic inconsistency (Haney, 1991), spurious flows and too much mixing. The fact that most, post-smoothing coastal depths in the present simulations were greater than 5 m would certainly affect bottom frictional drag and contribute to current inaccuracies. Of course, the greater number of nodes and elements in a refined grid would translate to larger computer memory requirements and longer run times. More recent versions of FVCOM than the one employed here are able to increase the time step size through the use of semi-implicit methods and this could at least partially offset the longer run times.

The sensitivity tests conducted with varying discharge rates and salinities for the Estero Basin underline the need for capturing the characteristics of rivers entering the region more accurately. Based on experience in the Broughton Archipelago, it is unlikely that the WSC will deploy new instrumentation on ungauged rivers without external financial support. Given the large costs involved and the financial limitations of the present project, this is probably not feasible. Thus, better methods for estimating discharges than those described in Section 2e need to be devised. Research has already begun on this and will be described in future manuscripts.

Though the present network of weather stations provides reasonable coverage of spatial variations in the winds and other atmospheric variables, compromises in site selection because of heavily forested terrain and steep-sided coastlines meant that some sites were partially sheltered from winds from particular directions. And although a credible technique was used for interpolating and extrapolating the observed winds to the entire model grid, missing and questionable data and smoothness assumptions inherent in this technique undoubtedly meant that the specified wind forcing had errors in some sub-regions. A high resolution atmospheric model that perhaps assimilates these observations would be an ideal solution to this problem but remains a challenge from both scientific and monetary perspectives.

Of course, the numerical model used in this application also has limitations in terms of its approximations to the underlying primitive equations and its parameterization of important processes such as mixing and bottom friction. No doubt, these limitations also contributed to the model inaccuracies seen in Figs 12, 13 and 14 and warrant further investigation that, it is hoped, will lead to improvements. Nevertheless, the foregoing model has provided a reasonably accurate first attempt at simulating flows in the Discovery Islands region that can be used for conducting particle tracking and dispersion studies. These could prove to be quite enlightening. Numerous Eulerian residual eddies, such as those displayed in Fig. 10, suggest the potential for retentive areas, but as demonstrated in Foreman et al. (1992) strong tidal flows may predominate, and the eddies could disappear in a Lagrangian sense. The results of these studies will be reported in future manuscripts.

Acknowledgements

We would like to thank Mainstream Canada, Marine Harvest Canada, Grieg Seafoods and Fisheries and Oceans Canada for partial financial support through project P-09-03-006 of the Aquaculture Coordinated Research and Development Program; Mainstream Canada for their temperature data from Barnes Bay and logistical support in deploying and maintaining weather stations in the region; the Program for Aquaculture Regulatory Research for partial funding; the Ecosystem Research Initiative of Fisheries and Oceans Canada for the funding to deploy mooring CM and Rick Thomson for permission to use its measurements; Tom Juhasz, Dave Spear and Lucius Perriault for deploying and recovering current meters; the Water Survey of Canada for providing river discharges; Trish Kimber for assistance with the figures; and two anonymous reviewers for constructive comments.