2.1 Model grid, parameterizations and parameter values, and particle tracking

The model used for this work evolved from the Finite Volume Community Ocean Model (FVCOM; Chen et al., 2003, 2006; Chen, 2022) application for the Discovery Islands developed by Foreman et al. (2012), herein referred to as F12. For the present work, the FVCOM code was upgraded from version 2.7 to 4.1 (Chen et al., 2013; Bianucci, 2024c). The grid was mostly the same as in F12 (horizontal resolution from ∼ 18 m in the narrowest channels to ∼ 1.2 km at the Strait of Georgia boundary; resolution calculated as the square root of twice the element area), except for a refinement in Bute and Toba inlets to improve the representation of the flow in those deep, narrow, and steep-sided fjords. The model bathymetry in these inlets was regenerated using a 10 m resolution digital elevation model (Kung, 2021), and the model grid excluded the steepest sections of the nearshore while also increasing the resolution, particularly in Bute Inlet (mean resolution of 182 m instead of 351 m as in F12). The resulting grid exhibited a total of 39 532 nodes and 72 518 elements (Fig. 11a). The model bathymetry was smoothed with a volume preserving technique that limits the ratio $\mathrm{\Delta }h/h\mathit{<}\mathrm{0.1}$ within each triangle in Bute and Toba inlets; everywhere else, the threshold was 0.3. The number of vertical terrain-following layers increased from 20 to 40, keeping higher resolution near the surface (surface layer represented 0.1 % of the total water column; the coarser layer at the bottom represented 7.5 % of the total depth).

Only a few parameterizations and parameter values changed in the new version of the FVCOM Discovery Islands model with respect to F12. The new setup used the k–ε vertical turbulence closure scheme (Rodi, 1987) and a background vertical diffusion and viscosity of 10⁻⁵ m² s⁻¹. The remaining parameterizations and parameters stayed the same and most are listed here for convenience. For instance, the Smagorinsky eddy parameterization was used for horizontal diffusivity with a coefficient C=0.02, and the horizontal viscosity was 10 times the diffusivity. Furthermore, the bottom roughness equation remained based on the General Ocean Turbulence Model (Burchard and Bolding, 2001) with a length scale of 10⁻³ m and a minimal value of $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ for the model bottom drag coefficient. The external and internal time steps were kept at 0.075 and 0.75 s (i.e., ISPLIT parameter set to 10). At the two open boundaries (one in the Johnstone Strait and the other in the Strait of Georgia, Fig. 1a), implicit radiation conditions were used for temperature and salinity (Blumberg and Kantha, 1985) and clamped conditions for surface elevation (Beardsley and Haidvogel, 1981); moreover, a sponge layer of 7.5 km wide was set with a damping coefficient of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{4}}$.

The Lagrangian particle-tracking module developed by Chen et al. (2006, 2013) was adapted at Fisheries and Oceans Canada to run offline, using hourly velocity outputs from FVCOM. This computationally efficient particle-tracking model is called PTrack (Losier, 2015) and simulates the particle trajectory until any of three termination conditions occur: advection outside the model domain, encountering land (“grounding”), or exceeding a user-imposed tracking limit (the latter was not implemented in this study). This Lagrangian model has been used for many applications related to tracking of pathogens, viruses, salmon post-smolts, etc. (e.g., Foreman et al., 2015; Quinn et al., 2022; DFO, 2022). Here, we followed the application of PTrack described in a study of the hydrodynamic connectivity between marine finfish aquaculture facilities in BC (DFO, 2022), using a time step of 10 s and outputting particle locations every hour.

2.2 Model forcing and initialization

The external forcing and initialization changed substantially from the previous versions of this model. First of all, the time period of simulation changed from April 2010 (F12) and April–October 2010 (Foreman et al., 2015) to May–June 2019, given the focus on understanding the stability of the subsurface temperature minimum found that year in Bute Inlet (Jackson et al., 2023).

Furthermore, the present setup benefited from new operational models to force the boundaries of the FVCOM domain. At the surface, hourly forcing of winds, pressure, and heat and water fluxes were provided by the High Resolution Deterministic Prediction System (HRDPS). The first release of HRDPS offered a resolution of 2.5 km (Milbrandt et al., 2016), which was suboptimal for a region with narrow inlets (∼ 1 km wide) such as the Discovery Islands area; therefore, for this study we used the experimental HRDPS version at 1 km horizontal resolution (Experimental HRDPS, 2022). The SalishSeaCast model (Soontiens et al., 2016; Soontiens and Allen, 2017), a three-dimensional physical–biological–chemical ocean model for the Strait of Georgia and Salish Sea, provided hourly temperature, salinity, and surface elevations for both the Johnstone Strait and Strait of Georgia open boundaries. Moreover, temperature and salinity fields from SalishSeaCast were used to initialize the FVCOM simulation, except in Bute Inlet, where observed profiles from 23 May 2019 allowed for the best possible initialization of the subsurface temperature minimum.

The 11 rivers included in the model (Fig. 1a) were implemented analogously to F12, except that the temperature and salinity at the discharge nodes were calculated using a mass conservation approach rather than specifying the actual value (i.e., the RIVER_TS_SETTING parameter was set to “calculated” instead of “specified”). As in the previous study, only 4 of the 11 rivers were gauged by the Water Survey of Canada during 2019 (Homathko, Campbell, Salmon, and Oyster rivers) and for the seven ungauged rivers, a watershed area-ratio approach was used to estimate their discharge. In F12, all ungauged rivers were estimated as the Homathko River discharge multiplied by their watershed area ratio (i.e., area of ungauged river divided by area of Homathko River). However, the new version of the model took advantage of a recent watershed characterization and classification (Giesbrecht et al., 2022) to select a representative donor gauge for each ungauged watershed. Specifically, each ungauged watershed was assigned a donor gauge from a nearby watershed of the same type and with similar climate and hydrology. Hence, while the Southgate and Toba rivers still used the Homathko River (all three having a glacierized mountain watershed type), the Stafford, Apple, Phillips, and Brem rivers based their discharges on the Wakeman River (snow mountain watershed type). The Powell River was estimated from the Campbell River, given that both are snow-mountain-type watersheds where discharge is controlled by dams (i.e., the underlying assumption being that both rivers were managed in the same way. The resulting river discharge dataset is available in Bianucci (2024b).

2.3 Model simulations

The model was run for ∼ 1 month, from 24 May to 27 June 2019. The choice of this period was determined by several factors. First of all, observations from summer 2019 in Bute Inlet showed evidence of the subsurface temperature minimum (and oxygen maximum) generated the previous winter during an Arctic outflow wind event (Jackson et al., 2023). Observations in Bute Inlet from 23 May (eight profiles) provided the temperature and salinity initial conditions in the fjord. Moreover, HRDPS-1 km outputs were available starting in 24 May 2019 (i.e., limiting any potential earlier start date). The total length of the simulation allowed for 5 d of spinup and 29 d for analysis; the latter is an appropriate averaging period to remove tides and calculate residual flows (Foreman et al., 1992). Longer simulations were not pursued partly because of the diffusive nature of FVCOM, which makes it challenging to reproduce a deep and narrow fjord without data assimilation. Nevertheless, the chosen month properly represented the summer conditions in the inlet, since the freshwater forcing was high from late May to late September, the wind conditions were stable (blowing mostly from the south until September), and the values at the open boundaries were similar throughout the summer (see Appendix B).

To represent idealized summer conditions in the absence of strong deep winter mixing the previous winter (e.g., by an Arctic outflow wind event), a sensitivity experiment was performed by removing the temperature minimum feature in Bute Inlet from the initial conditions. This experiment represents an extreme scenario, given that strong Arctic outflow wind events are common in winter in this region (more than two events per year on average; Jackson et al., 2023), such that some degree of subsurface cooling is usually present (e.g., MacNeill, 1974; Pickard, 1961). All other initial conditions and forcings (e.g., atmospheric and open boundaries) remained unchanged, given that the winter deep-mixing event would only affect summer conditions in the fjord (e.g., summer open boundary conditions in the Strait of Georgia and Johnstone Strait would not be affected by the outflow winter event in Bute Inlet). It is unclear how the winter event might have affected the summer river discharge; we kept this forcing unchanged to focus on the role of the initial conditions, acknowledging that this assumption is a source of uncertainty. In this sensitivity simulation, initial temperature and salinity profiles in Bute Inlet were kept constant below the main pycnocline (Fig. A1); the constant values corresponded to the coolest and saltiest observations in the deepest third of the water column.

The same particle-tracking experiment was run with velocity outputs from each simulation (referred to as “baseline” and “sensitivity” depending on whether they were initialized with the observed profiles in Bute Inlet or not, respectively). Virtual particles were released at the start of both simulations, using the same initial locations in both cases. Namely, particles were released within the location of the observed cold subsurface feature, i.e., at every grid node in Bute Inlet where temperature was ≤ 8 °C in the baseline initial conditions. Particles were tracked during the whole length of the simulations as they were advected by the 3D flow fields; if particles reached the coastline or seafloor at any given time, they became grounded and were removed from the particle count. Time series of the non-grounded particles remaining in Bute Inlet were calculated. We had the ability to let particles trapped in the bottom bounce back into the water column if the bottom vertical velocity was upwards (i.e., include particle resuspension); however, results did not change significantly and are not shown.

2.4 Available observations and metrics for model evaluation

Observed vertical profiles were available for model evaluation from bottle and conductivity–temperature–depth (CTD) measurements (Bianucci, 2024a). In total, 114 profiles were available for the period 24 May–27 June (see locations in Fig. 1b), representing 20 348 matches between modelled and observed temperatures (20 231 matches for salinity). Sea surface elevations were available at the Campbell River tidal gauge (triangle in Fig. 1b). Unfortunately, no observed velocity profiles were available in Bute Inlet to evaluate the modelled currents.

To evaluate the performance of the model, potential temperature (θ) and salinity fields were compared against observed values. For this purpose, model outputs were selected from the grid node closest to each observation and linearly interpolated in the vertical dimension and time to create model–observation pairs for each available in situ sample. With all the available pairs, we calculated several metrics frequently used to quantify model–observations misfit (both for the full model domain and for the top 100 m only). The model bias determines the mean deviation between modelled and observed values, while the root mean square error (RMSE) measures the deviation in a least-squares sense. These two metrics retain the units of the analysed variable. Two useful non-dimensional metrics are the model efficiency (ME) or skill and the Willmott skill score (Willmott, 1981). The former relates the deviation between model and observations to the variability in the observations, while the latter is an indication of the model error divided by the range of the observations. Both skill metrics can range from zero to one, with one indicating perfect agreement between observations and model. The latter is also the case for the square of the Pearson's correlation coefficient (R²), which quantifies the correlation between modelled and observed data. A succinct but thorough description of these metrics can be found in Lehmann et al. (2009) and Liu et al. (2009).

4.1 Baseline simulation with observed initial conditions

Temperature and along-inlet velocity were averaged over the last 29 d of the simulation. This averaging effectively filtered out the tides (to focus on residual flows), while also removing the first 5 d of simulation to allow for spinup. A transect plot of mean along-inlet velocities through Bute Inlet showed a multi-layered structure of the velocity field in most of the fjord (Fig. 4a). The surface layer flowed outwards of the fjord, with a return flow underneath down to approximately the depth of the outside sill (∼ 300 m), following a typical estuarine circulation. However, the return flow had a clear vertical structure, with velocities close to zero at the depth of the minimum averaged temperature (Fig. 4a and b; a dotted horizontal line at 50 m highlights the co-location of the near-zero averaged velocities and the mean temperature minimum). Below the depth of the outside sill, the mean, slow flow was towards the mouth of the inlet, with a narrow and weak inflow layer near the seafloor. At around 70 km from the head of the inlet, mean velocities showed the effect of tidal mixing over the outer sill (Fig. 4a) given the strong tidal currents in the Discovery Islands region (Foreman et al., 2012, 2015).

Figure 4Mean along-inlet transects throughout Bute Inlet for (a, b) baseline and (d, e) sensitivity simulations. Variables shown are (a, d) mean along inlet velocity and (b, e) mean potential temperature. Velocities are positive (red) towards the mouth of the inlet and negative (blue) towards the head. Averages over the last 29 d of the simulation removed tidal effects. The dotted horizontal line at 50 m highlights the location of the mean temperature minimum in all panels. (c) Map of Bute Inlet transect, colour-coded by the distance from the head of the inlet.

To further explore the vertical structure of potential temperature and along-inlet velocity, profiles were plotted every 10 km in the middle of the inlet (from 20 to 50 km away from the head of the inlet; Fig. 5a to d). Four or five layers are evident in the mean along-inlet velocity profiles when considering the return flow between ∼ 5 and ∼ 300 m depth to be divided in two where the velocities approach zero, i.e., above and below the depth of the temperature minimum feature. The temperature minima in each profile are found between 44 and 46 m deep, with mean velocities there ranging between −0.01 and −0.7 cm s⁻¹ (negative values imply flow towards the head of the fjord).

Figure 5Vertical profiles of mean along-inlet velocity (coloured red/blue) and potential temperature (grey) for the (a–d) baseline and (e–h) sensitivity simulations, at four locations in the inlet (from left to right: 20, 30, 40, and 50 km away from the head; see Fig. 4c). Velocities are positive (red) towards the mouth of the inlet and negative (blue) towards the head. Velocity values for the outflowing surface layer are given as red numbers on the top-right of each panels (in cm s⁻¹).

Download

4.2 Sensitivity simulation without temperature minimum feature

A sensitivity simulation with homogeneous conditions under the main pycnocline (as described in Sect. 2.3) showed a single return flow from underneath the surface outflow to the depth of the outside sill (Figs. 4d and 5e to h) instead of the two return layers found at those depths on the baseline simulation. The surface, outward-flowing layer was almost identical in both simulations (Figs. 4 and 5), and the bottom outward and weak inward near-seafloor flows were quite similar (e.g., comparable magnitude, vertical distribution/shape, and depth range). Ignoring the weak bottom mean inflow between 40 and 60 km, the vertical mean flow structure in this simulation could be described as three-layered (Fig. 5e to h).

The cold subsurface feature can be identified as a new water mass when comparing temperature–salinity diagrams for both the baseline and sensitivity simulations (Fig. 6a; only model values at the time and location of the observations were plotted to simplify the figure). The strong winter mixing event not only led to subsurface cooling and reoxygenation but also affected salinity (see profiles in Fig. 3c, d) and led to density and stratification changes when comparing both simulations (Fig. 6c, d). In particular, stratification decreased below the outward-flowing layer between ∼ 5 and 50 m in the baseline (negative values in Fig. 6d) and density increased in that same region (positive values in Fig. 6c). These changes were not uniform in the horizontal but stronger near the head of the inlet (seen in Fig. 6c and d and highlighted by the Δρ profiles at 20 and 40 km in Fig. 6b); given the overall horizontal density gradient near the surface (lower density near the head due to the freshwater inputs), the subsurface water mass led to a reduced horizontal density difference between the head and the mouth of Bute Inlet. The latter in turn led to an even slower mean along-inlet velocity underneath the estuarine surface outflow in the baseline simulation.

Figure 6(a) Temperature–salinity diagram for the baseline (red) and sensitivity (black) simulations in Bute Inlet. Model results are shown at the time and location of the observations (12 and 26 June 2019; location in Fig. 4g); for reference, four isopycnals were labelled according to their σ_θ (kg m⁻³). (b) Profiles of mean density difference between both simulations (Δρ) at 20 and 40 km away from the head of the inlet, shown in the top 100 m of the water column. Bottom panels show along-inlet transects of the difference in (c) mean density and (d) mean Brunt–Väisälä frequency (N²) between baseline and sensitivity experiments (Δ=29 d average baseline minus 29 d average sensitivity; negative values indicate that the baseline is less dense/stratified than the sensitivity simulation).

Download

4.3 High retention of subsurface particles

The slow velocities below the surface outflow seen in both the baseline and sensitivity simulations (Figs. 4a, d and 5) also led to high retention in Bute Inlet in both particle-tracking experiments (Fig. 7). The time series of the percent of particles moving inside the fjord (i.e., not grounded) showed that more than 96 % of the moving particles stayed within the inlet (north of 50.45° N) by the end of the simulations. Particularly, the retention was higher in the baseline simulation (thicker line in Fig. 7), with more than 97 % of the particles remaining in the inlet. Furthermore, the median depth of the particles in the baseline simulation stayed within 45 and 54 m (comparable to the initial median depth of 54 m), contrasted with the deeper median depth of the sensitivity simulation (> 100 m after 5 d, thinner line in Fig. 7). The weaker stratification below the main pycnocline led to a larger range of particle dispersion and more grounding in the latter simulation.