Seasonal Variability of the Ocean Circulation in Queen Charlotte Strait, British Columbia

ABSTRACT

Queen Charlotte Strait (QCST) is a large marine area separating northern Vancouver Island from the British Columbia mainland. Although it is an important waterway connecting Queen Charlotte Sound with Johnstone Strait and eventually the Strait of Georgia, complex dynamics in QCST, including tides, surface winds, river runoffs, and coastal currents, are yet understudied. In this numerical study, a high-resolution (up to approximately 10 m in the horizontal) model was developed for QCST and adjacent areas to investigate the circulation in the strait, using the unstructured grid, Finite-Volume, primitive equation Community Ocean Model (FVCOM). Two operational larger-scale models force the FVCOM model: the Coastal Ice Ocean Prediction System for the West coast (CIOPS-W) at the open ocean (one-way nesting), and the High-Resolution Deterministic Prediction System (HRDPS) at surface boundaries. Prevailing winds over QCST are southeasterly in winter and northwesterly in summer and the model demonstrated the capacity to reproduce the ocean circulation regimes in both seasons of 2019. Complex empirical orthogonal function (CEOF) and correlation analyses determined that surface winds play a dominant role on the variability of residual flow in QCST in the winter, while river runoffs strongly influence the variability during the summer. The analyses highlight that QCST is a complex estuarine system with significant seasonal variability rather than just a simple conduit between Queen Charlotte Sound and the Strait of Georgia.

RÉSUMÉ

[Traduit par la rédaction] Le détroit de la Reine-Charlotte est une vaste zone marine qui sépare le nord de l’île de Vancouver du continent de la Colombie-Britannique. Bien qu’il s’agisse d’une importante voie navigable reliant le détroit de la Reine-Charlotte au détroit de Johnstone et éventuellement au détroit de Géorgie, la dynamique complexe du détroit de la Reine-Charlotte, y compris les marées, les vents de surface, les écoulements fluviaux et les courants côtiers, est encore peu étudiée. Dans cette étude numérique, un modèle à haute résolution (jusqu’à ∼10 m à l’horizontale) a été élaboré pour le détroit de la Reine-Charlotte et les zones adjacentes afin d’étudier la circulation dans le détroit, en utilisant le modèle océanique communautaire à grille non structurée, à volume fini et à équation primitive. Deux modèles opérationnels à plus grande échelle forcent le modèle océanique communautaire à grille non structurée, à volume fini et à équation primitive: le Coastal Ice Ocean Prediction System for the West Coast (CIOPS-W) au niveau de l’océan ouvert (imbrication unidirectionnelle), et le High-Resolution Deterministic Prediction System (HRDPS) aux frontières de surface. Les vents dominants sur le détroit de la Reine-Charlotte sont du sud-est en hiver et du nord-ouest en été et le modèle a démontré sa capacité à reproduire les régimes de circulation océanique dans les deux saisons de 2019. Les analyses de fonction orthogonale empirique complexe et de corrélation ont déterminé que les vents de surface jouent un rôle dominant sur la variabilité du débit résiduel dans le détroit de la Reine-Charlotte en hiver, tandis que les écoulements fluviaux influencent fortement la variabilité pendant l’été. Les analyses soulignent que le détroit de la Reine-Charlotte est un système estuarien complexe avec une variabilité saisonnière importante plutôt qu’un simple conduit entre le détroit de la Reine-Charlotte et le détroit de Géorgie.

KEYWORDS:

• Queen Charlotte Strait
• FVCOM
• coastal modelling
• tidal current
• residual flow
• wind-driven current
• river runoffs
• complex empirical orthogonal function

1 Introduction

Queen Charlotte Strait (QCST) in the British Columbia (BC) central coast is a large marine area separating northern Vancouver Island from mainland BC. It is an important waterway connecting Queen Charlotte Sound (QCS) with Johnstone Strait and eventually the Strait of Georgia; moreover, it connects QCS with many coastal inlets through the Broughton Archipelago (Fig. 1(a)). This coastal corridor of BC is a region with significant biodiversity and the host to a highly productive ecosystem. The complex coastal waterways support a vibrant fish farming industry, i.e. numerous salmon farms are located in QCST (Haggarty et al., 2003), and provide subsistence and traditional ways of life for many First Nations communities. The route through Johnstone Strait and QCST is used by Fraser River sockeye salmon smolts when they leave the Strait of Georgia for the Pacific Ocean, as well as by adult sockeye salmon when they return from the ocean (Quinn & Terhart, 1987).

Fig. 1 Location of Queen Charlotte Strait (QCST) and general seasonal surface circulation in the region (Figure adapted from Borstad et al., 2011). Winds and surface currents off Vancouver Island are generally (a) southeasterly in the winter and (b) northwesterly in the summer. Ocean shading denotes water depths as indicated on the colour scale. The maps highlight geographic locations mentioned in the main text.

QCST is characterised by complex coastlines due to mainland inlets and islands, and it receives freshwater discharges from many rivers. It is a relatively broad pass compared with Johnstone Strait, with a considerably greater cross-sectional area. QCST is approximately 90 km long and its main basin widens from around 13 km to 26 km (Thomson, 1981). The central strait reaches depths of 400 m (Crawford et al., 2007) and connects with Cook Trough in QCS (Fig. 1(b)). Goletas Channel (between Hope Island and Vancouver Island) is a deep channel (approximately 350 m) along the northwest of Vancouver Island, restricted by a 140 m sill (Nahwitti Bar) at the western end of the channel (Fig. 1(a)).

Ocean currents in the BC central coast are mainly forced by tides, surface winds, and freshwater discharges (Thomson, 1981). The tide floods from south QCS into QCST with a large tidal range of 5.5 m as observed in Port Hardy (Crawford et al., 2007). The surface currents in QCST can generate a short, steep chop more than a metre or two high when winds rise above 10 m s⁻¹ or 20 knots (Thomson, 1981). Seasonal patterns in the wind-driven currents have been found off BC’s central coast, in association with the large-scale atmospheric pressure systems offshore and consequently the regional wind forcing (Crawford et al., 2007). Southeasterly winds during autumn and winter drive shelf surface currents eastward through the Scott Islands towards the coast (Crawford et al., 2007; Freeland et al., 1984). These southeasterly winter winds relax in the spring and become weaker and more frequent northwesterly winds in summer (Borstad et al., 2011). These northwesterly winds drive summer upwelling at the eastern QCS and wind-driven surface currents exit QCS at Cook Bank (Crawford et al., 1985). This feature has been supplemented by analysis of satellite imagery data sets (Crawford et al., 2007; Han & Chen, 2022) and tracks of Argos surface drifters (Crawford et al., 1999). Figure 1 shows a schematic of the general seasonal circulation in the region.

As part of the regional wind system, prevailing winds in QCST are also southeasterly in winter and northwesterly in summer. Coverage and resolution from satellite data is poor for QCST. Present knowledge of the flow structure in QCST region is based on limited historical measurements. Waters in QCST are well mixed or weakly stratified at wintertime, while waters are stratified in the summer due to increased runoffs and thermal heating at surface (Stucchi et al., 2002). Many large watersheds in the region are snow/glacier fed, leading to larger discharges in summer (Giesbrecht et al., 2022). Among them, river runoffs from Rivers Inlet and Smith Inlet on the northwest (see Fig. 1(b)) are the most significant source of fresh water (Table 1). Waters along the north shore of QCST are colder on average, because of discharges of rivers fed by snow/ice melt, and possibly because of upwelling and tidal mixing (Crawford et al., 2007). Based on observed water properties (e.g. the salinity distribution near the surface), the large volumes of fresh water contributed by rivers tends to move at the surface along the north shore into QCS (Crawford et al., 1999; Scagel, 1961). Estuarine circulation in QCST drives surface water out of the strait (Stucchi et al., 2002; Thomson, 1981). The oceanic water apparently flows into the Strait along the deep channels in the central part of the Strait, as well as along Vancouver Island (Scagel, 1961). Therefore, the northwest end of QCST, connecting with QCS, represents a complicated transition between shelf and coastal dynamics.

Table 1. Rivers included in the model and their mean discharges during the modelling periods.

Rivers	Name	Location	Watershed Type	Mean Discharge (m^3s^−1)
				Dec 26, 2018–Jan 31, 2019 (historical Jan)	Jul 6–Aug 12, 2019 (historical Jul)
1	Klinaklini	East	Glacierized	76 (80)^a	696 (732)^a
2	Glendale	East	Snow	11	16
3	Kakweiken	East	Snow	23	34
4	Ahta	East	Snow	17	26
5	Kingcome	East	Snow/Rain	73	109
6	Wakeman	East	Snow/Rain	56 (49)^b	84 (99)^b
7	Tsitika	Vancouver Is.	Rain	43 (33)^c	7 (10)^c
8	Nimpkish	Vancouver Is.	Rain	78	13
9	Wannock	Northwest	Snow/Glacierized	186 (171)^d	411 (587)^d
10	Bella Coola	Central BC	Snow/Glacierized	32 (36)^e	231 (245)^e
11	Dean	Central BC	Snow	4 (5)^f	28 (24)^f
12	Nekite	Northwest	Snow/Glacierized	93	206
13	Keogh	Vancouver Is.	Rain	13 (12)^g	2 (1)^g
Total				574	2302

Note: Watershed type follows characterisation by Giesbrecht et al. (2022). For the seven gauged rivers, the observed historical mean discharge for January and July is provided in parenthesis; superscripts indicate the years of observations available for each river gauge: (a) 1977–2020, (b) 2006–2020, (c) 1975–2020, (d) 1928–2020, (e) 2012–2019, (f) 1960–2021, (g) 2016–2020.

The circulation of large areas of the BC northern and central coast has been studied using numerical modelling methods. However, QCST was often set as an open boundary in those models (e.g. Crawford et al., 1999; Cummins & Oey, 1997; Foreman & Thomson, 1997). QSCT was fully covered in Foreman et al. (2004) using a 2D unstructured model to study tidal dissipation in the waters around Vancouver Island. Foreman et al. (2006) and Foreman, Czajko, et al. (2009) studied the estuarine and tidal circulation in the Broughton Archipelago (to the southeast of QCST) using the Finite Volume, primitive equation Community Ocean Model (FVCOM (Chen et al., 2006; Chen, Beardsley, Cowles, et al., 2013; Chen, Beardsley, Luettich, et al., 2013)). Given the focus on the Broughton Archipelago, this model did not cover QCS (to solve deep water exchanges) nor Rivers Inlet and Smith Inlet (to include the significant river runoffs into QCST). Khangaonkar et al. (2017) also created a FVCOM model with a domain completely encompassing Vancouver Island. However, since their area of interest was Puget Sound (southern end of the Salish Sea, within Washington State), errors in elevations at Port Hardy within QCST were relatively high.

A better understanding of the synthesis of the tidal flow, wind-driven circulation, freshwater discharges, and deep-water exchanges in the QCST area would assist the aquaculture industry and its regulators, particularly as BC moves towards an Area-Based Aquaculture Management (ABAM) approach (Fisheries and Oceans Canada, 2022). For this purpose, a high-resolution numerical ocean model for QCST was developed using FVCOM. Here, the model was validated for winter and summer of 2019, using in situ observations of currents, temperature, and salinity. Analyses of model results focus on seasonal flow regimes in western QCST and investigate the roles of surface winds and river runoffs on the residual flow variability.

The paper is organised as follows: Section 2 describes the model development, Section 3 reports on model performance, and Section 4 describes the mean and tidal flows and examines the variations of residual flow and winds. Finally, Section 5 discusses the results and conclusions of the study.

2 FVCOM model development for QCST

In this numerical study, a high resolution FVCOM model was developed with a focus on QCST. The model domain extends from QCS on the northwest and to Johnstone Strait on the southeast (Fig. 2). The triangular unstructured grid has 123,598 nodes with 221,232 elements. The horizontal grid resolution (square root of element area) varies from approximately 9.6 m for small islands and channels to about 1.9 km in the open shelf to the west. In the vertical, a sigma-coordinate system with 20 levels has higher resolution near both the surface and the bottom and relatively lower resolution at mid-column; more explicitly, inter-layer spacing ranges from 0.03% of the total water depth at surface and 1.5% at bottom to 15% at mid-column. For a typical water depth of 100 m, there are 11 layers within the upper 20 m approximately, allowing for the resolution of the pycnocline and mixing layer in the summer. The grid was generated with Gmsh (version 2.5), developed by Geuzaine and Remacle (2009).

Fig. 2 Model mesh and bathymetry (depths greater than 500 m in dark red). The location of the 13 rivers included in the model is shown, as well as geographical features mentioned in the text. The greyed area along the western boundary marks the nesting zone. Mooring sites in central Queen Charlotte Strait (LBA-1a and LBA-1b) and Goletas Channel (LBB) are marked with letters A, C, and B, respectively.

The area around Goletas Channel required some of the highest horizontal resolution to solve the narrow, deep, and steep-sided channel along with its entrance sill. The model grid also resolved the detailed north shore of QCST by including the Belize and Seymour inlets (Fig. 2), where the tidal currents can reach 8 m s⁻¹ or 16 knots at the narrows connecting to QCST (Dallimore & Jmieff, 2010).

The 25 m resolution Canadian Hydrographic Service (CHS) high water mark lines were used to generate the coastlines of the model domain. The lines provide the level reached by sea water at high tide as used for CHS nautical charts. The British Columbia 3 arc-second Bathymetric Digital Elevation Model (Lemieux & Luquire, 2021) was used to generate the model bathymetry. To reduce the hydrostatic inconsistency associated with sigma-coordinate hydrodynamic modelling, the model bathymetry was smoothed with a volume preserving technique that limits the ratio Δh/h < 0.3 within each triangle (e.g. Hannah & Wright, 1995), where Δh is the difference between maximum and minimum depths at the three nodes of the triangle and h is the average depth. Minimum water depths were set to 5 m in order to eliminate wetting and drying in the model grid.

The FVCOM model version 4.1 (Chen et al., 2006; Chen, Beardsley, Cowles, et al., 2013; Chen, Beardsley, Luettich, et al., 2013) employs a numerical integration with a second-order accurate finite volume flux discrete scheme, where internal and external modes are ‘split’ and integrated over distinct time steps. The external time step was set to 0.075 s to preserve stability over the simulation. The ratio of time steps between the internal and external modes was set to 10 in this application. The Smagorinsky eddy parameterisation was used for horizontal diffusivity with coefficient C = 0.1. The model was closed mathematically using the Generalised Ocean Turbulent Model (GOTM (Burchard et al., 1999; Umlauf & Burchard, 2003)) with the k-ϵ scheme for vertical mixing and a vertical Prandtl number of 1.0 (i.e. viscosity and diffusion used the same vertical mixing coefficient). The background vertical diffusion and viscosity were set to 10⁻⁴ m²s⁻¹. The bottom roughness equation was based on the GOTM as well, with a length scale set to 0.002 m and a minimal value of 0.005 for the model bottom drag coefficient.

The external forcing for the coastal ocean model was facilitated by Canadian regional operational forecast systems. At the surface, the High Resolution Deterministic Prediction System (HRDPS (e.g. Milbrandt et al., 2016)) provided hourly surface winds, pressure, and heat and water fluxes at a ∼2.5-km horizontal resolution. Heat fluxes were computed with standard bulk formulae, using surface water temperatures from FVCOM itself and air temperature, air pressure, humidity, cloud cover, and short and long wave radiation from HRDPS fields. The Coastal Ice Ocean Prediction System for the West Coast of Canada (CIOPS-W (e.g. Lu et al., 2017; Paquin et al., 2021)), a Northeast Pacific Ocean model with a horizontal resolution of 1/36° (approximately 2–2.5 km), provided information at the two open boundaries of the FVCOM domain. One of the boundaries is at the east of the domain in Johnstone Strait (Fig. 2), where hourly surface elevations and daily temperature and salinity profiles provided lateral boundary conditions. These profiles were assumed to be uniform over the cross-section of the eastern open boundary due to the coarse resolution of CIOPS-W over the narrow Johnstone Strait. The second boundary lays at the open shelf, from the north of QCS at Price Island to Brooks Peninsula on western Vancouver Island. A 15-node-wide nesting zone was implemented along this boundary to better represent shelf and coastal currents as well as open ocean waters. Inside the nesting zone, FVCOM nodes matched the model grid from CIOPS-W (greyed area in Fig. 2) and FVCOM model variables (temperature, salinity, and currents) were one-way nested from CIOPS-W hourly fields with a timescale of an hour. A hyperbolic tangent weight function of distance was applied within the nesting zone, varying from one at the outer boundary (i.e. fully diagnostic) to zero at the inner boundary (i.e. fully prognostic). Hourly surface elevations from CIOPS-W were forced at the open boundary nodes, such that FVCOM tidal forcing is limited to the same eight main constituents used in CIOPS-W (M2, K1, N2, S2, K2, O1, P1, Q1), plus the associated overtides. The sum of the eight constituent amplitudes is approximately 85% of the maximum tidal range as calculated at tidal gauges around Vancouver Island (M. G. G. Foreman, personal communication, November 1, 2022).

The model included 13 rivers (Fig. 2 and Table 1); daily discharge observations from the Canadian Hydrological Database (wateroffice.ec.gc.ca) were available for seven of those rivers (Klinaklini, Wakeman, Tsitika, Wannock, Bella Coola, Dean, and Keogh). The approach discussed by Foreman, Czajko, et al. (2009) was used for five of the six ungauged rivers (Glendale, Kakweiken, Ahta, Kingcome, and Nimpkish), as these rivers were included in their model for the Broughton Archipelago. For the remaining ungauged river (Nekite), half of the Wannock River discharge was used in the FVCOM model, the same ratio used in CIOPS-W (Paquin et al., 2021). Furthermore, FVCOM requires time series of river temperature and salinity as part of its freshwater inputs. Given the lack of such riverine observations, river temperatures were set to the sea surface temperature values from the CIOPS-W outputs closest to the river mouth. For river salinity time series, two approaches were used depending on the representation of the rivers in the FVCOM grid: (1) those rivers with mouths explicitly represented in the grid (rivers 1–8 in Fig. 2) were assigned zero salinity, and (2) rivers flowing directly at the coastline nodes (rivers 9–13 in Fig. 2) were assigned surface salinities from CIOPS-W minus 6 (the subtracted value was selected through calibration). Overall, all but three of the 13 rivers are fed by snow and/or glacier melt, leading to peak discharges in summer (Giesbrecht et al., 2022); the three exceptions are small rain-fed rivers on Vancouver Island. The largest rivers in the region are the Klinaklini at the east end of the model domain (via Knight Inlet) and Wannock in QCS (via Rivers Inlet; see locations in Fig. 2 and discharges in Table 1).

Three model simulations were performed for this study (Table 2). Model results from simulations Winter-V and Summer-V&A (January 19 to February 28, 2019, and July 1 to August 11, 2019, respectively) were used for model validation and evaluation, based on the availability of the Acoustic Doppler Current Profiler (ADCP) and Conductivity, Temperature, and Depth (CTD) observational data described in Section 3. Model results from simulations Winter-A (December 21, 2018 to January 31, 2019) and Summer-V&A were used for analysis, because those time periods corresponded to prevailing southeasterly winds in winter (downwelling-favourable; Fig. 3(a)) and northwesterly winds in summer (upwelling-favourable; Fig. 3(b)). Climatological winds (35-year reanalysis data, January 1, 1979–October 2, 2014) from the North American Regional Reanalysis (NARR (Mesinger et al., 2006)) by the National Centers for Environmental Prediction (NCEP) were compared with time-averaged HRDPS winds (red and black arrows, respectively, in Fig. 3). Winter 2019 winds were typical, blowing in similar directions and magnitudes as in the climatology. Summer 2019 wind patterns suggest that the North Pacific High pressure system was weaker and further south than in the climatology, which is consistent with previous findings (e.g. Amaya et al., 2020). Nevertheless, the largest differences were found offshore, such that winds over QCST and southern QCS in summer 2019 were relatively similar to climatological conditions. Furthermore, freshwater discharges in gauged rivers were also close to historical mean values in winter and summer 2019 (Table 1), such that results from the winter and summer 2019 simulations represent relatively normal/typical conditions in QCST.

Fig. 3 Temporal average of HRDPS wind fields shown in black arrows for winter (a, December 26, 2018–January 31, 2019) and summer (b, July 6–August 11, 2019), respectively. Climatological winds from NARR for the same two periods shown in red vectors (climatologies based on years 1981–2010).

Table 2. Description of model simulations.

Model Runs	Spin-up Period	Examined Period	Description
Winter-V	Jan 19–23, 2019	Jan 24–Feb 28, 2019	Winter simulation for validation
Winter-A	Dec 21–25, 2018	Dec 26, 2018–Jan 31, 2019	Winter simulation for analysis
Summer-V&A	Jul 1–5, 2019	Jul 6–Aug 11, 2019	Summer simulation for validation and analysis

All model simulations were integrated for 41-day periods, starting from initial temperature and salinity fields from CIOPS-W (velocities were set to zero) and allowing for a spin-up period of 5 days. Model results of the last 36 days were examined and presented in this study, using hourly outputs of water levels, currents, water temperature, and salinity for model evaluation and analysis.

3 Model validation

The model performance was evaluated using measurements available in QCST during winter and summer 2019, including tidal gauge data from CHS, and ADCP and CTD observations by Fisheries and Oceans Canada (DFO) obtained from www.waterproperties.ca (publicly available, but registration required; also, publicly available at www. cioos.ca without registration). The following subsections describe the comparison of the model against each of these three datasets.

a Port Hardy Tidal Gauge

Modelled water levels and derived tidal constituents at Port Hardy were compared against derived tidal gauge constituents, obtained by analysing observed surface elevations for the same time periods as the model simulations (Table 3). The tidal analysis used Codiga’s Unified Tidal Analysis and Prediction programme UTide (Codiga, 2011). This programme is a suite of modern tidal analysis and prediction algorithms implemented in Matlab® and built on the work of previous software packages (e.g. Foreman, 1977, 1978; Foreman, Cherniawsky, et al., 2009; Leffler & Jay, 2009; Pawlowicz et al., 2002).

Table 3. Tidal constituent comparison between model and observational tidal gauge data at Port Hardy.

Model Runs	Tidal Constituent	Amplitude (m)			Phase (°)			Complex Distance (cm)
		OBS	MOD	DIFF	OBS	MOD	DIFF
Winter-V	M2	1.32	1.26	−0.06	253.49	252.68	−0.82	6.51
	K1	0.53	0.51	−0.01	267.87	265.98	−1.89	2.21
	S2	0.46	0.41	−0.04	297.17	296.35	−0.82	4.53
	O1	0.29	0.29	0.00	231.89	232.84	0.95	0.48
	N2	0.33	0.26	−0.07	229.32	230.07	0.75	6.88
	Q1	0.05	0.04	−0.01	226.70	226.64	−0.06	0.96
Summer-V&A	M2	1.32	1.25	−0.07	253.86	252.90	−0.96	7.43
	K1	0.60	0.60	0.01	263.80	261.14	−2.66	2.87
	S2	0.34	0.35	0.01	294.77	296.27	1.51	1.59
	O1	0.30	0.29	−0.01	233.49	234.16	0.67	0.73
	N2	0.29	0.26	−0.03	239.90	229.71	−10.19	5.70
	Q1	0.05	0.04	−0.01	236.51	225.13	−11.38	1.13

Major tidal constituents in this region are M2, K1, S2, O1, N2, and Q1, listed in descending order of tidal amplitudes. In the comparison between model and observational tidal gauge results, the model reproduced well both tidal amplitudes and phases. The error of each tidal constituent was computed as the discrepancy in the complex plane between the model results and tidal gauge derived data (complex distance column in Table 3, calculated as in Foreman et al. (2012)). The largest complex distances for all six tidal constituents corresponded to M2 and N2 tides (6.51 and 6.88 cm, respectively, in winter and 7.43 and 5.70 cm, respectively, in summer). Models with complex distance values within approximately 5% of the observed amplitude are generally considered to have acceptable accuracy, which is the case for all harmonics except S2 in winter (11%) and N2 and Q1 in both seasons (21–28%). The two major constituents (M2 and K1) were properly represented within the expected 5% threshold. To improve the modelled tides (particularly S2, N2, and Q1), boundary elevation conditions (currently from CIOPS-W) could be tuned, particularly in the narrow Johnstone Strait (not well resolved by CIOPS-W).

b ADCP

Two moorings with ADCP instruments near the surface and bottom were deployed by DFO on January 21, 2019, in central Queen Charlotte Strait (LBA-1a and LBA-1b) and Goletas Channel (LBB; see locations in Fig. 2). The upward-looking ADCPs had a bin size of 8 m in the vertical; after post-processing, ocean currents at a 30-minute interval from January to August 2019 were available. LBA-1a was deployed at 50° 53.996′N, 127° 38.410′ W in 448 m of water; only the near-bottom current metre data were used because of flotation issues on the upper ADCP. This mooring was hit and relocated 110 m shallower to 50° 56.121′N, 127° 41.173′W on February 26, 2019, and remained in the new location until August 10, 2019; it was renamed as LBA-1b for this period. LBB was deployed at 50° 50.519′N, 127° 42.998′ W in 396 m of water until August 11, 2019, with two upward-looking ADCP units at approximately 332 m (near-bottom) and 36 m. Therefore, current metre measurements from bottom to subsurface (roughly 315 m to 27 m) and subsurface to near-surface (approximately 35 m to 5 m) were available for both winter (model simulation Winter-V) and summer (model simulation Summer-V&A) periods.

Model generated eastward and northward components (U and V) of near-bottom ocean currents at the LBA-1a and LBA-1b sites were compared with ADCP data (Fig. 4). Observed near-bottom (approximately 422 m) currents in winter were dominated by tidal currents and displayed large amplitudes in both U and V components (up to about 0.5 m s⁻¹; red lines in Fig. 4(a,b)). Model generated ocean currents agreed with the winter ADCP data with a root-mean-square error (RMSE) of less than 0.15 m s⁻¹ in both U and V components. During the 2019 summer simulation, the model also performed well in central QCST at LBA-1b, with RMSE values of 0.12 m s⁻¹ (U) and 0.09 m s⁻¹ (V) at near surface (around 25 m), and 0.10 m s⁻¹ (U) and 0.07 m s⁻¹ (V) near-bottom (approximately 313 m). The differences between modelled and observed currents are likely partly from residual flow and/or due to the lack of minor tidal harmonics (and associated tidal range) in the FVCOM forcing (since CIOPS-W only includes eight main constituents).

Fig. 4 Current component comparison between the ADCP observations (red curves) and model results (blue curves) for (a) U and (b) V near-bottom at LBA-1a (Fig. 2) in winter, (c) U and (d) V subsurface at LBA-1b in summer, and (e) U and (f) V near-bottom at LBA-1b (Fig. 2) in summer.

Near surface currents measured by the top ADCP at the LBB site demonstrated strong non-tidal variations in both winter (Fig. 5(a,b)) and summer (Fig. 5(e,f)). Eastward current components (U) are much greater than northward flow (V) due to the narrow geometry of Goletas Channel. Overall model generated near surface currents are comparable with ADCP data, with RMSE values of 0.16–0.18 m s⁻¹ for U components (Fig. 5(a,e)) and 0.08 m s⁻¹ for V components (Fig. 5(b,f)). The relatively large discrepancy between modelled and observed surface U components was driven mainly by winter differences between January 27 and February 3. This discrepancy was found at surface and decreased rapidly with depth (e.g. RMSE values for U decreased from 0.15 m s⁻¹ at 27 m to 0.09 m s⁻¹ at 50 m in LBB’s near-bottom ADCP, Table 4), suggesting that the HRDPS wind forcing at ∼2.5 km resolution is likely not high enough to drive properly the surface flow in the narrow Goletas Channel (approximately 2 km wide). Wind observations available later in the year at Bull Harbour in Goletas Channel confirmed that, at times, HRDPS generated stronger west or northwest winds than those observed by the weather station (not shown). For the near-bottom currents, the model-observations agreement improved for all components of the flow, with RMSE between 0.04 and 0.08 m s⁻¹ (Fig. 5(c,d) for winter; g and h for summer). Both summer and winter simulations produced a good representation of the tidal and non-tidal features of the flow at LBB.

Fig. 5 Current component comparison between the ADCP observations (red curves) at the LBB site (Fig. 2) and model results (blue curves) for (a) U and (b) V near surface and (c) U and (d) V near-bottom in winter, and (e) U and (f) V near surface and (g) U and (h) V near-bottom in summer. Near surface and near-bottom observations correspond to the top and near-bottom ADCPs, respectively.

Table 4. Root mean square errors (RMSE) of modelled currents (winter and summer validation simulations, Table 2) vs. near-bottom ADCP observations at different depths as well as for the depth averaged flow (LBA-1a not included because only the bottom bin was usable). The top bin (∼25–27 m) represents the near-surface flow as sampled by the near-bottom ADCP.

Periods	Winter-V		Summer-V&A
Near-bottom ADCP	LBB		LBA-1b		LBB
Depths	RMSE (U, m s^−1)	RMSE (V, m s^−1)	RMSE (U, m s^−1)	RMSE (V, m s^−1)	RMSE (U, m s^−1)	RMSE (V, m s^−1)
Top bin: LBA-1b:25 m; LBB:27 m	0.15	0.07	0.12	0.09	0.15	0.07
50 m	0.09	0.06	0.14	0.14	0.11	0.06
100 m	0.10	0.05	0.08	0.10	0.08	0.05
150 m	0.07	0.04	0.08	0.08	0.06	0.05
200 m	0.06	0.04	0.08	0.10	0.06	0.05
250 m	0.05	0.04	0.08	0.12	0.06	0.05
300 m	0.07	0.05	0.09	0.08	0.05	0.05
Bottom bin: LBA-1b:313 m; LBB:315 m	0.08	0.05	0.10	0.07	0.05	0.04
Averaged RMSE of all depths	0.08	0.05	0.09	0.10	0.07	0.05
RMSE of Depth Averaged Flow	0.03	0.02	0.04	0.04	0.03	0.02

Overall, the high resolution 3-D model demonstrated the ability to simulate both winter and summer flow regimes at all depths, performing particularly well at representing the depth-averaged flow (summary of velocity RMSE scores shown in Table 4, both for different locations in through the water column and for the depth-averaged flow).

c CTD

Moored CTD data were collected by DFO in central QCST (LBA-1a and LBA-1b) and Goletas Channel (LBB) (locations A, C, and B, respectively, Fig. 2). The CTD observations every 30 min covered the same January to August periods as the ADCP moorings. At LBA-1a, continuous measurements were available at three depths: subsurface (40 m), middle depth (148 m), and near-bottom (436 m). The LBA-1b site only provided near-bottom (327 m) CTD data, while CTD sampling at LBB in Goletas Channel occurred at subsurface (36 m) and near-bottom (332 m).

The model reproduced the observed temperature and salinity at LBA-1a, LBA-1b, and LBB well (Fig. 6). The linear regressions showed slope values close to 1 (i.e. unbiased) for both modelled temperature and salinity. The RMSE varied between 0.27 and 0.54°C for temperature and 0.17 and 0.31 for salinity (values shown in Fig. 6).

Fig. 6 Modelled (MOD) (a) temperature and (b) salinity compared against CTD data (OBS) at subsurface (36 m) and near-bottom (332 m) at the LBB site in winter and summer, subsurface (40 m), middle depth (148 m), and near-bottom (436 m) at the LBA-1a site in winter, near-bottom (327 m) at the LBA-1b site in summer. The regression slope and RMSE are listed in the same order as in the legend. CTD mooring sites are marked in Fig. 2.

Temperature and salinity CTD profiles inside QCST during the winter and summer simulation periods were available from a few cruises by Fisheries and Oceans Canada (DFO) and compared against model results (Fig. 7). There were 18 CTD profiles collected at 15 sites during the winter period and 7 profiles/sites available during the summer period (blue crosses and red dots, respectively, in Fig. 7(c)). In winter, observed water temperature and salinity ranged from 7.4°C and 33.2 in the deep waters (more than 400 m) to 9.3°C and 30.6 at surface. In summer, observed values ranged from 6.6°C and 33.8 at depth (more than 400 m) to above 12.4°C and 31.5 at surface within QCST. Generally, the ocean exhibited stronger stratification in the summer. The model captured the winter temperature and salinity well, with RMSE of 0.57°C for water temperature and 0.67 for salinity fields. In the summer, temperature RMSE increased to 0.77°C, mainly caused by the overestimation of sea surface temperature at some sites. However, the modelled salinity during summer improved on the winter metrics, with an RMSE value of only 0.33.

Fig. 7 (a) Temperature and (b) salinity scatter plots using the CTD profiles (OBS) and model results (MOD); blue and red colours indicate winter and summer data, respectively. Panels show RMSE and regression slope values. (c) Locations of the CTD sites in the study area during the winter (blue crosses) and summer (red dots) simulation periods.

In conclusion, the high resolution 3-D model was able to represent appropriately the observed temperature and salinity in the two seasons, both at the moorings and the profiling stations.

4 Ocean circulation in QCST

This section describes in detail the modelled circulation in QCST, particularly investigating the variability of the ocean currents as influenced by tides, winds, and freshwater discharges during both winter and summer.

a Mean Seasonal Flow Regimes

Surface and near-bottom flows from simulations Winter-A and Summer-V&A were time-averaged to produce mean seasonal patterns (Fig. 8). To remove any tidal bias, both simulations were averaged over 29 days (December 30, 2018, to January 28, 2019, and July 10 to August 8, 2019). In the winter, when the temporally averaged surface wind fields blew to the northwest (Fig. 3(a)), the surface water mean flow was northwestward, joined by the currents coming from the west of Vancouver Island around the Cook Bank area (Fig. 8(a)). Meanwhile, the deep currents headed into QCST to form a typical estuarine structure in the vertical (Fig. 8(c)). The deep near-bottom flow was intensified in the narrow deepening in the central QCST; such an intrusion of near-bottom water along the deep channels in the central part of QCST was also noted by Scagel (1961). Depth in Goletas Channel is also large but the shallower Cook Bank area at the western entrance constrains deep-water intrusions from QSC and/or from around northern Vancouver Island.

Fig. 8 Mean surface and near-bottom (1 m above seabed) flows based on the (a, c) Winter-A and (b, d) Summer-V&A simulations. The averaging period was 29 days to ensure the removal of tidal bias (December 30 to January 28, 2019, and July 10 to August 8, 2019). Near-bottom flow only shown where bathymetry was greater than 50 m.

In summer, the average surface flow was weaker than in winter (Fig. 8(b)) due to the reversal of the temporally averaged winds, which blow into QCST in summer (Fig. 3(b)). Surface waters moved northwestward mainly along the north shore, as previously described by Crawford et al. (1999), driven by the large summer freshets of mainland rivers (e.g. the Klinaklini River within our model domain, see Table 1). Furthermore, a time-mean southeastward flow along the south shore of QCST formed an apparent counterclockwise pattern at surface (Fig. 8(b)), consistent with observations by Scagel (1961) and Crawford et al. (2007). At depth, the time-mean deep flow in QCST was remarkably similar to winter in location, direction and strength (Fig. 8(d)), suggesting that (1) changes in the surface seasonal drivers (winds and runoff) may compensate each other at depth and/or (2) external, non-seasonal drivers are at play (e.g. from the shelf boundary conditions). Despite the seasonal differences at surface and similarities at depth, the estuarine flow structure in QCST was maintained.

b Tidal Flow

Tidal analysis was applied to simulated surface current speeds in both winter and summer using UTide (Codiga, 2011). We focus the discussion of the model results on the region between 50.76 and 51.2°N and between 128.4 and 127.3°W; this area in western QCST had the highest horizontal resolution (a modelling choice driven by the location of aquaculture interests in the region). As the M2 and K1 tidal constituents have the largest amplitudes in this region, their tidal current ellipses derived from model results over QSCT were illustrated at the surface layer for both seasons (Fig. 9). For M2, the strongest tidal flow was found at the narrows connecting the Belize and Seymour inlets with the north of QCST (Fig. 9(a,b)). Modelled tidal M2 currents were enhanced by local topographic features around the islands in central QCST and around the northern Vancouver Island coast. There, tidal flow from the broad shelf squeezes into Goletas Channel over the sill (Nahwitti Bar, Fig. 2), leading to increased tidal kinetic energy due to the topographic influence of the local narrow channel. For the second largest tidal component, K1, surface tidal flow is enhanced only at the entrance to Goletas Channel (Fig. 9(c,d)). Surface tidal ellipses for both harmonics showed no evidence of major seasonal differences between winter and summer, except some differences in central QCST (Fig. 9).

Fig. 9 Tidal ellipses (M2 and K1) of surface flow in the winter (a, c) and summer (b, d) 2019 simulations.

The near-bottom M2 and K1 tidal current ellipses (1 m above seabed) were also derived from model results, but only for areas with water depths greater than 50 m (Fig. 10). The major semi-axis magnitudes for the tidal currents are larger at surface and weaker near-bottom (Fig. 9(a,b) and Fig. 10(a,b)). For M2, large near-bottom tidal flow was found at the deepest part of central QCST and in Goletas Channel (Fig. 10(a,b)). This type of localised baroclinic (internal) tides was not found for K1 tidal flow in the QCST area (Fig. 10(c,d)). As in the case of the surface, generally there was no large seasonal differences in the near-bottom tidal flows (Fig. 10).

Fig. 10 Tidal ellipses (M2 and K1) of flow at 1 m above seabed (where water depths are greater than 50 m) in the winter (a, c) and summer (b, d) 2019 simulations.

c Non-tidal Flow: CEOF Analysis on Winds and Residual Currents

Complex empirical orthogonal function (CEOF) analyses on both the winds and currents allow a further examination of the role played by seasonal winds on the circulation variability in QCST. Vectors were transformed into a complex scalar during the analysis process. Winds from HRDPS (10 m above sea surface) were used for this analysis (same wind fields used to force the FVCOM model). For the currents, tidal flows (e.g. Fig. 9) were removed to focus on the residual component of the circulation. The time series of both wind and current vectors were low-pass filtered with a timescale of 48 h to focus on synoptic variability. The winter and summer analysis periods were defined as in the Winter-A and Summer-V&A simulations: December 26, 2018 to January 31, 2019 (when the surface is dominated by southeasterly, downwelling-favourable winds, Fig. 3(a)), and July 6 to August 11, 2019 (when the surface winds are northwesterly and upwelling-favourable, Fig. 3(b)). CEOF analyses were conducted on a structured grid with a horizontal resolution of 2 km (both HRDPS winds and FVCOM current fields were interpolated to this new structured grid). The Belize and Seymour inlets in the north of QCST were removed from the analysis because of the coarse resolution (2 km) of the structured grid.

For the winds, the percentages of the variance accounted for by the first two CEOF modes were 96.47% and 3.02% in winter and 89.72% and 9.35% in summer. Within QCST, the 1st CEOF mode represented northwesterly winds in the winter and southeasterly winds in the summer (Fig. 11(a,b)). Furthermore, the amplitude of the 1st CEOF mode in summer was about half of the one in winter. The 2nd mode, already minor compared with the 1st CEOF mode, illustrated an anticyclonic pattern in the winter (centred at the middle of the entrance of QCST; Fig. 11(c)), while showing a diverging pattern in the summer (separating at the entrance of QCST; Fig. 11(d)).

Fig. 11 Amplitudes (units: m s⁻¹; background colours) and vectors of the first two CEOF modes based on low-pass filtered (48 h cut-off) wind fields at 10 m height from HRDPS model results during (a and c) winter (December 26, 2018, to January 31, 2019), and (b and d) summer (July 6 to August 11, 2019).

For the surface currents, the first two winter CEOF’s explained 87.94% and 8.52% of the variance, explaining together 96.46% of the variance of the winter surface residual flow (Fig. 12(a,c)). In summer, the first two CEOF’s explained 76.61% and 18.16% of the variance (together accounting for 94.77%) of the surface residual flow in QCST (Fig. 12(b,d)).

Fig. 12 Amplitudes (units: m s⁻¹; background colours) and vectors of the first two CEOF modes based on low-pass filtered (48 h cut-off) model reproduced surface currents during (a and c) winter (December 26, 2018, to January 31, 2019), and (b and d) summer (July 6 to August 11, 2019).

The 1st CEOF modes of surface residual flow in winter and summer showed uniform directions throughout the study region (Fig. 12(a,b); vector axes come out of the CEOF analysis, but vector directions were chosen to match those of the mean seasonal flows). In winter, the surface residual current direction patterns were similar to the 1st CEOF mode of winds (Fig. 11(a)), which is consistent with how the mean surface flow responds to the prevalent northwesterly winds (Fig. 8(b)) and suggests that QCST is part of the coastal current system driven by the large-scale winds (Fig. 3). In summer, mode 1 current directions were opposite and not fully aligned to those of mode 1 winds (Fig. 11(b)), implying that other mechanisms must primarily drive the surface flow variability. In terms of the amplitude of the variance, the 1st CEOF modes of surface residual flow showed strong variance in the central QCST and Goletas Channel in both winter and summer (Fig. 12(a,b)). In summer, strong variance was also found along the north shore of QCST. This feature could be related to the runoff from the mainland inlets along the north shore and at the east end of QCST given the high river discharge in summer (Table 1) and the low salinity zone found there in the summer simulation (not shown). The latter low salinity zone has also been observed in previous studies (e.g. Crawford et al., 2007; Scagel, 1961).

The directions of the 2nd CEOF mode of surface residual flow in winter (Fig. 12(c)) were similar to the 2nd CEOF mode of winds during the winter (Fig. 11(c)), with an anticyclonic pattern centred at the middle of the entrance of QCST. Variations were strong in the middle of QCST, Goletas Channel, and the Cook Bank (Fig. 1(b)). As shown in Fig. 12(c), Cook Bank is the place where the currents coming from the west of Vancouver Island meet the outflow variance from Goletas Channel and the north of Hope Island (Fig. 1(a)).

In the summer, directions of the 2nd CEOF mode of surface residual flow (Fig. 12(d)) were quite different compared with the 2nd CEOF mode of winds (Fig. 11(d)), suggesting that different dynamics affect the secondary variance of the surface residual flow within QCST in summer. For example, a cross-strait residual flow pattern within QCST suggests the presence of upwelling-driven variability along the north shore, and downwelling-driven variability at central and the south shore of QCST. Furthermore, a strong variance zone along the north shore of QCST could be enhanced by the runoff and associated low surface salinity/density there.

d Surface Residual Flow Variability in QCST

The corresponding principal component (PC) time series (units: non-dimensional) of CEOF modes for winds and surface residual flow were strongly correlated in winter (Fig. 13). The correlation coefficients were 0.87 for CEOF mode 1 and 0.74 for CEOF mode 2, both significant at the 99% confidence level (degrees of freedom of the low pass filtered time series is 35). These correlations suggest that the surface residual flow anomalies are mainly driven by winds in winter; i.e. the strong wind variance along the northwest-southeast axis drives the circulation associated with the 1st mode of circulation (Fig. 11(a) and Fig. 12(a)), while the anticyclonic gyre over the QCST is linked to the mode 2 responses (Fig. 11(c) and Fig. 12(c)). Nevertheless, although the directions of the surface residual currents were consistent with the CEOF maps of the winds, the variance strength in space was modified by the topographic features in QCST (Fig. 12(a,c)).

Fig. 13 The corresponding principal component (PC) time series (units: non-dimensional) for (a) 1st CEOF mode and (b) 2nd CEOF mode of low-pass filtered (with a cut-off timescale of 48 h) surface winds (blue lines) and surface residual flow (red lines) during the winter period.

In contrast, no significant correlation between PC time series of winds and surface residual flow existed in the summer (time series not shown), with correlation coefficient values of −0.13 (CEOF mode 1) and −0.27 (CEOF mode 2). These low correlation values suggest that winds are not the main factor affecting the surface residual flow patterns shown in Fig. 12(b,d). The latter is further supported by the different spatial patterns of vectors from wind CEOF modes (Fig. 11(b,d)) and surface residual flow CEOF modes (Fig. 12(b,d)).

River runoff, and consequently, ocean stratification variance/water density anomalies, could play an important role on surface residual flow variance, especially in the summertime when large discharges occur (Table 1). The largest river flows belong to glacierized watersheds, whose discharges peak in summer due to glacier melt: Klinaklini in Knight Inlet, Wannock in Rivers Inlets, and Nekite in Smith Inlet (Fig. 2). Nevertheless, these runoffs must cover a long path to reach QCST, such that tidal and wind-driven processes modify the freshwaters along their way. Therefore, time series of volume discharges at the river mouth were found to be a poor indicator of the effect of runoff on surface residual flow variance (not shown).

To investigate the role of rivers on the variance of summer surface residual flow variability, modelled surface salinity was used as the indicator. Correlation coefficients between surface salinity (low-pass filtered with a cut-off timescale of 48 h) and PC time series of CEOF modes 1 and 2 of surface residual flow (i.e. PC corresponding to Fig. 12(b,d)) were calculated at each model node and mapped (Fig. 14). The area of analysis was expanded (compared with the area of CEOF analysis) to include more regions with sources of fresher water (e.g. outflows from Rivers Inlet to the north). Negative correlations suggest a physical meaning here, since lower surface salinity likely implies larger river outflow (neglecting contributions from precipitation), which leads to stronger surface residual flow in the directions shown in Fig. 12(b,d) (higher salinity leads to stronger residual flows in opposite directions). Surface residual flow could also drive salinity changes through upwelling/downwelling and lead to negative correlations; however, the low correlation between winds and surface residual flow in summer suggests that wind-driven upwelling/downwelling is not a main driver of the residual circulation. Significant negative correlation was found to the north in the region connecting with Rivers Inlet for CEOF mode 1, while for CEOF mode 2 significant negative correlation was along the north shore of QCST, a region influenced by the rivers to the east.

Fig. 14 Correlation coefficients between surface salinity (low-pass filtered with a cut-off timescale of 48 h) and PC time series corresponding to CEOF (a) mode 1 and (b) mode 2 of surface residual flow (CEOF modes shown in Fig. 12(b,d)). Dashed boxes mark the area of the CEOF analysis of currents.

To identify the potential dynamics and from where the low salinity waters (negative correlation) come from, a series of correlation maps were created with time lags (salinity leading) ranging from 2.5 to zero days for CEOF mode 1 (i.e. strong flow along the north shore, Fig. 12(b)), and from 1.5 to zero days for CEOF mode 2 (i.e. the crossing strait mode, Fig. 12(d)) of surface residual currents. The significant correlation for CEOF mode 1 was tracked to Rivers Inlet (Wannock River) with a time lag of 2.5 days (Fig. 15). The significant correlation for CEOF mode 2 originated from rivers located in eastern QCST via the pathway along the north shore (Fig. 16) with a shorter time lag (∼1 d). These results strongly suggest that the negative correlations respond to riverine influence, since lower surface salinity was tracked to freshwater sources, leading to higher variance in surface currents. Lagrangian tracer simulations could help to support these findings in future studies. The role of upwelling/downwelling other than wind-driven (e.g. wave-driven) remains to be studied.

Fig. 15 Correlation coefficients between surface salinity (low-pass filtered with a cut-off timescale of 48 h) and PC of CEOF mode 1 of surface residual flow (Fig. 12(b)), with different time lags (salinity leading 2.5 days to 0 days). Dashed boxes mark the area of the CEOF analysis of currents.

Fig. 16 Correlation coefficients between surface salinity (low-pass filtered with a cut-off timescale of 48 h) and PC of CEOF mode 2 of surface residual flow (Fig. 12(d)), with different time lags (salinity leading 1.5 days to 0 days). Dashed boxes mark the area of the CEOF analysis of currents.

e Bottom Residual Flow Variability

A similar CEOF analysis was performed on the deep flow by examining the near-bottom residual flow variability at 1 m above the seabed. The details of the CEOF analysis were analogous to that of the surface residual currents, e.g. tides removal, 48-hour low-pass filter, 2 km resolution structured grid, and periods of analysis corresponding to winter and summer.

The first two winter CEOF modes explained 71.15% and 19.98% of the variance, respectively, explaining together 91.13% of the variance for the bottom residual flow in the winter simulation (Fig. 17(a,c)). In the summer, the first two CEOF’s only explained 56.75% and 25.73% of the variance (a total of 82.48%) of the bottom residual flow in QCST (Fig. 17(b,d)).

Fig. 17 Amplitudes (units: m s−¹; background colours) and vectors of the first two CEOF modes based on low-pass filtered (48 h cut-off) model reproduced near-bottom currents (1 m above seabed) during (a and c) winter (December 26, 2018, to January 31, 2019), and (b and d) summer (July 6 to August 11, 2019).

The amplitude and vectors of the 1st CEOF mode of near-bottom residual flow in winter illustrated a clear path for the deep-water exchange between QCS (via Cook Trough, Fig. 1) and QCST via the central deep channel in QCST (Fig. 17(a)). However, the amplitude of the variance was much less than surface residual flow (Fig. 12(a)). Furthermore, while this deep flow path was visible in the 1st CEOF mode for the summer (Fig. 17(b)), the amplitude of the variance was even smaller than in winter (less than 1 cm s⁻¹).

In the temporal domain, the PC times series of the 1st CEOF mode of the low-pass filtered near-bottom residual flow were significantly correlated with surface wind variance in both winter and summer (Fig. 18(a,b)). In winter, positive correlation coefficients were found between southeasterly winds (CEOF mode 1, Fig. 11(a)) and both surface/near-bottom residual flow pattern (CEOF mode 1, Fig. 12(a) and Fig. 17(a); surface correlation in Fig. 13).

Fig. 18 The corresponding PC time series (units: non-dimensional) for 1st CEOF mode of low-pass filtered (with a cut-off timescale of 48 h) surface winds (blue lines) and near-bottom residual flow (1 m above seabed, red lines) during (a) the winter simulation period and (b) the summer simulation period.

In summer, a significant positive correlation was found between northwesterly winds (CEOF mode 1, Fig. 11(b)) and near-bottom residual flow pattern (CEOF mode 1, Fig. 17(b)); i.e. when the winds blow into QCST, the near-bottom residual currents flow northwestward and leave the Strait.

Although the 2nd CEOF mode accounted for 19.98% of the variance of the near-bottom residual flow in winter and 25.73% in summer, the spatial distribution was very patchy (Fig. 17(c,d)). These spots of large variance of the near-bottom residual currents were likely related to localised dynamics such as topographic effects and baroclinic instabilities. They were very localised and constrained by the local topographic features (e.g. bottom slope). For example, relatively large amplitude is found near the sill (Nahwitti Bar, Fig. 2) in Goletas Channel (Fig. 17(d)).

5 Discussion and conclusions

Queen Charlotte Strait is an important waterway with significant biodiversity and the host to a highly productive ecosystem. It is a meeting place for the relatively weakly stratified waters at Johnstone Strait that move seaward and the more highly stratified shelf waters that move inland from Queen Charlotte Sound (Thomson, 1981). Because of the large-scale atmospheric pressure systems offshore, and their resulting regional wind forcing, there is a strong seasonal nature in the wind-driven current patterns in QCST. Furthermore, the role of river runoffs is significant in summer because many large watersheds in the region are snow/glacier fed with summer freshets.

Using a numerical ocean model, this study analysed the seasonal circulation in QCST, investigating mean, tidal and residual flows. Model simulations were examined for representative periods corresponding to prevailing southeasterly winds in winter (December 26, 2018, to January 31, 2019) and northwesterly winds in summer (July 6 to August 11, 2019). Wind and freshwater forcing in both study periods of 2019 were mostly typical, except for the summer offshore structure of the winds that responded to a weaker North Pacific High pressure system. Nevertheless, given that the analyses focused on local winds in the south QCS and QCST regions, the results are likely generalisable to any regular winter and summer season in QCST.

In winter, the temporally averaged surface flow was northwestward (same direction as the mean winds), joined by the currents coming from the west of Vancouver Island around the Cook Bank area (Fig. 8(a)). Meanwhile, an intrusion of near-bottom water along the deep channels in the central part of QCST formed a typical estuarine structure in the vertical (Fig. 8(c)). In summer, the average surface flow was weaker than in winter and mostly against the surface winds, with the strongest flow along the north shore (Fig. 8(b)); however, the mean deep flows remained like in winter (Fig. 8(d)). Therefore, the estuarine flow structure in QCST was maintained in summer by large seasonal runoffs (particularly from rivers on the mainland of BC with large summer freshets) and in winter by a combination of smaller runoff and wind forcing. Future work should explore the role of external (e.g. open boundaries) and internal dynamics leading to similar mean near-bottom flows in QCST in winter and summer.

Complex empirical orthogonal function (CEOF) and correlation analyses revealed the influence from seasonal winds and freshwater discharges. In winter, variability of the surface residual flow was mainly along-strait (87.94% of the variance, CEOF mode 1), which was similar to the dominant direction of the winter winds (96.47% of the variance, CEOF mode 1). The large portion of the variance explained by mode 1, in addition to the significant correlation between the associated PC times series of CEOF mode 1 (r = 0.87), strongly suggest that the variability of the surface residual flow is mainly driven by winds in the winter. The 2nd CEOF modes for wind and surface residual flow fields in winter also exhibited a similar spatial pattern (an anticyclonic gyre centred at the middle of the entrance of QCST); while the variance explained by the mode 2 was lower than for mode 1 (3.02% for winds and 8.52% for surface residual flows), the PC time series were still highly correlated (r = 0.74). In total, 96.46% of the variance of surface residual flow (CEOF modes 1 and 2) can be explained by winds in the winter, highlighting the main role that atmospheric forcing plays in the variability of the circulation in QCST in this season.

For the summer period, the analysis suggested that winds are not a significant driver of surface residual flow variability in QCST; in contrast, river runoff plays the most important role. The increased relative importance of the 2nd CEOF mode of surface residual flow (76.61% and 18.16% of variance explained by mode 1 and 2, respectively) suggests that the stronger ocean stratification brings more complex patterns to the total variance of the residual flow. Furthermore, the spatial patterns of vectors from wind CEOF mode 1 and from surface residual flow CEOF mode 1 were reversed in direction, with weakly correlated PC time series (correlation coefficient of −0.13). In contrast, correlations between the PC time series of surface residual currents with low-pass filtered surface salinity (as a proxy for riverine influence) were strong. In particular, the strong residual flow anomalies along the north shore in CEOF mode 1 (Fig. 12(b)) were tracked to Wannock River discharges (via Rivers Inlet) on the northwest, highlighting the large spatial extent of the influence of this river. The crossing strait residual flow anomalies from the north shore in CEOF mode 2 (Fig. 12(d)) are mainly trackable to the rivers to the east of QCST, likely from the inlets ending at the Broughton Archipelago and potentially from sources beyond Johnstone Strait. River runoffs, and consequently, ocean stratification variations and water density anomalies, play an important role on surface residual flow variance during the summer.

The CEOF analysis of the near-bottom residual flow variability at 1 m above seabed showed a deep-water exchange route between QCS and QCST via Cook Trough and QCST’s central deep channel in both seasons (mode 1, Fig. 17(a,b)); the variance of this deep flow represented less than 20% of the variance at surface (Fig. 12(a,b)). The significant correlations between the mode 1 near-bottom PC times series and those for winds (Fig. 18) suggest that winds play an important role in the deep-water exchange between QCST and QCS. Nevertheless, the decreased relative importance of the 1st modes compared with the 2nds suggest a more complex structure of the total variance of near-bottom residual flows (Fig. 17). Freshwater inputs still contribute to the estuarine circulation even in winter, likely driving part of the near-bottom residual flow variability.

In summary, main drivers of the residual circulation differ in winter and summer, according to the described CEOF analyses. In winter, the significant positive correlations between winds and both surface and near-bottom residual flows suggest that winter winds drive residual flow variations at all depths; however, the flow directions of the responses to winds are different at surface and near-bottom. When the winds blow towards the northwest, the surface residual currents flow northwestward and leave the Strait, while the residual deep circulation flows into the Strait (mainly through the central QCST). The latter suggests that the near-bottom wind-driven currents come from deep QCS to compensate for the surface outflow; this flow might be partly affected by the coastal current system driven by large-scale winds and freshwater/estuarine forcing. In contrast, the summer analysis suggests that surface flow variance is mainly driven by freshwater runoff; fluctuations of near-bottom currents likely result from a combination of the mass balance driven by the estuarine circulation as well as the coastal current system driven by large-scale winds.

In conclusion, this numerical study shows that QCST is a complex estuarine system rather than just a simple conduit between QCS and the Strait of Georgia, with variability driven by surface winds in winter and affected by river runoffs in summer. While some of the seasonal mean patterns in the general QCS/QCST area were previously described in other studies (e.g. Crawford et al., 1999; Freeland et al., 1984; Scagel, 1961; Thomson, 1981), QCST was rarely the region of focus. A synthesis of tidal flows, wind-driven circulation, freshwater discharges, deep water exchanges, and associated variability was lacking in this area. However, it is important to improve our understanding of the physical oceanography of QCST, since this region holds many fisheries and aquaculture operations and it is actively used for transport as well as for recreational and traditional purposes. In particular, given the importance of the region as a migration route of Fraser River sockeye salmon, future efforts could aim at using the model to understand interannual variability as well as climate change effects expected in the region, including biogeochemical impacts.

Acknowledgements

This work was partially funded by the Program for Aquaculture Regulatory Research (PARR), which is designed to increase our scientific knowledge and inform regulatory decision making and policy development. We gratefully acknowledge Dr. Gary Borstad, who provided copyright permission to use their Fig. 1 from Borstad et al. (2011). The authors thank Glenn Cooper, Mike Foreman, Maxim Krassovski, Jennifer Jackson, Hayley Dosser, and Pramod Thupaki for their contributions and discussions at different stages of this work. In particular, Mike Foreman provided a thoughtful review of the first draft of the manuscript. We also thank Michael Dunphy, Maxim Krassovski, and Pramod Thupaki for FVCOM code and tools, and David Spear for ADCP data collection and his expertise on data interpretation. Two anonymous reviewers provided valuable feedback.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by funding from Fisheries and Oceans Canada (DFO).