Impact of advection loss due to wind and estuarine circulation on the timing of the spring phytoplankton bloom in a fjord

Abstract

A coupled biophysical model is used to explore the physical controls involved in the timing of the spring phytoplankton bloom in fjords. Observations from Rivers Inlet, British Columbia, are used to force and evaluate the model. It is found that the interannual variation in timing is due primarily to variations in retention, in particular, to variations in horizontal advection out of the fjord. The two dominant processes are (i) strong outflow winds rapidly advecting the surface layer and thus the phytoplankton population out of the fjord and (ii) losses due to high river flux increasing the estuarine circulation. Both processes delay the timing of spring bloom. Smaller effects on the interannual variation are due to increased wind mixing which deepens the mixing layer and reduces light to phytoplankton, and increased river flow which increases the stratification and decreases the mixing layer depth. Observed interannual variations in cloudiness were small. Strong outflow winds are common in winter along the British Columbia coast, but generally cease after the spring wind transition. Thus, observed interdecadal variations in the spring transition date probably imply strong variations in the timing of spring phytoplankton blooms in British Columbia fjords.

biophysical model, fjords, interannual variation, spring bloom

Introduction

In temperate regions, phytoplankton follow a strong seasonal cycle which is typically characterized by a prominent spring phytoplankton bloom. The spring bloom occurs when increases in biomass of phytoplankton due to processes such as growth and accumulation exceed decreases in biomass of phytoplankton due to a variety of processes (Longhurst, 1995). Physical mechanisms directly or indirectly contribute to both types of processes. For example, growth rates of phytoplankton are determined by light and nutrient limitations and are mediated by water temperature (Parsons et al., 1984; Reynolds, 2006). The light level that a population of phytoplankton receives depends on the incoming light levels determined by time of year, latitude, cloud amount and albedo, the extinction of the light in the water, and by the mixing layer depth, that depth to which the phytoplankton are actively mixed by the surface turbulence (Sverdrup, 1953; Collins et al., 2009).

Sverdrup (1953) proposed a model in which physical processes affect the timing of the spring bloom by varying phytoplankton growth rates through varying the mixing layer depth. While this classical critical depth theory has been well developed and validated with observations in many situations including estuaries and in coastal upwelling areas (Mann and Lazier, 2006), the influence of the mixing layer depth may be of less importance in other marine environments (Lucas et al., 1998; Huisman et al., 1999). A number of studies have found variations in the strength of the loss terms to be dominant. A model of the spring bloom in shallow estuaries found that leakage of phytoplankton from the shallow surface layer due to sinking or turbulent diffusion in permanently stratified estuaries is important and must be taken into consideration when modelling the spring bloom (Lucas et al., 1998). Decreases in the strength of grazing have been shown to contribute to the timing of the start of the spring bloom in the Atlantic (Behrenfeld et al., 2013). Other factors affecting biomass loss are advection (Collins et al., 2009) and viral lysis (Welschmeyer and Lorenzen, 1985).

In other systems, variability in accumulation rather than growth or loss rates is the dominant bloom initiation process. For example, in northerly Norway and West Spitsbergen, the onset of the spring phytoplankton bloom is regulated by the accumulation of biomass due to the entrainment of phytoplankton from resting spores (Eilertsen et al., 1995). As the transition from resting spore to active phytoplankton is determined by daylength, in high latitude fjords, the spring phytoplankton bloom can occur at approximately the same time each year (Eilertsen et al., 1989). In the open ocean, deep mixing can change the effective daylength leading to interannual variation in the transition to active phytoplankton and thus in the timing of the bloom (Mignot et al., 2014).

In temperate British Columbia fjords and in our Rivers Inlet example in particular, we expect the interannual variation in the timing of the spring bloom to be dependent on mixing layer depth, advective loss, and incoming radiation. Winter concentrations of zooplankton are similar between years and concentrations do not increase significantly until after the spring phytoplankton bloom occurs (Tommasi et al., 2013b).

British Columbia is well known for its complex coastline, interwoven with fjords and inlets. Fjords are associated with high latitudes and glacial activity. They have steep sides and deep, elongated channels with a sill that traps deeper water. Fjord circulation is generally characterized by moderate to large river discharge and weak to moderate tidal forcing due to the large depths, resulting in a strongly stratified estuary (Pedersen, 1978; Valle-Levinson, 2010). A shallow surface layer transports freshwater out towards the inlet mouth while deeper and often nutrient-rich waters enter the inlet over the sill near the mouth of the fjord. Winds are generally funnelled by the steep fjord sides to be either inflow winds (towards the head) or outflow winds (towards the mouth). The strong stratification in the fjord allows the winds to generate strong flows above the halocline. Sustained winds large enough to reverse the surface layer flow frequently occur in Knight Inlet, a British Columbia fjord (Baker and Pond, 1995). Both estuarine circulation (Kaarvedt and Svendsen, 1990) and wind events (Winter et al., 1975; Furuya et al., 1993) have been shown to advect and dilute plankton populations in fjords and bays.

Rivers Inlet is a temperate coastal fjord, located on the central coast of British Columbia (Figure 1). Historically, it hosted one of the largest sockeye salmon runs in Canada until stocks collapsed in the 1990s (McKinnell et al., 2001). The inlet is uniformly deep (200–300 m), is 40 km long and 3 km wide, and has a relatively deep and complex sill region at its mouth (depth 110 m), opening into Queen Charlotte Sound. The major source of freshwater is from Owikeno Lake which drains into the Wannock River which, in turn, drains into the head of the inlet.

Figure 1.

Map of Rivers Inlet on the West Coast of British Columbia, including sample station DF02 (diamond), the Florence Daily site (circle), the Wannock River Mouth (end of line), the Laska and Ethel meteorological stations (black x), and Port Hardy airport and Cathedral Point Automated Weather Station (filled circles).

There is considerable interannual variability in the timing and duration of the spring bloom in Rivers Inlet (Tommasi et al., 2013a). In other regions, it has been suggested that the timing of the spring bloom may be influential to the survival of various higher trophic levels through a match/mismatch process (Cushing, 1974, 1975). For example, in the Strait of Georgia, early spring blooms can cause nauplii of Neocalanus plumchrus, which migrate to the surface in spring, to miss the spring bloom causing them to rely solely on post-bloom phytoplankton which are less abundant and less nutritious (El-Sabaawi et al., 2009). In Rivers Inlet, a late spring bloom resulted in a smaller, delayed zooplankton biomass peak and a change in the zooplankton community composition (Tommasi et al., 2013a). Zooplankton are the primary food source of juvenile sockeye salmon (Foerster, 1968; Pauley et al., 1989; Ajmani, 2012), and thus changes in the timing, biomass, and composition of these lower trophic levels may be an important factor in sockeye early marine survival, and ultimately fluctuations in the sockeye salmon stock biomass.

We use a coupled biophysical model (Collins et al., 2009) to determine which physical forcings most strongly influence the timing of the spring phytoplankton bloom in Rivers Inlet. A quasi one-dimensional vertical mixing model is used to model the physical properties of Rivers Inlet. The physical model is a K-Profile Parameterization (KPP) non-local boundary layer model (Large et al., 1994) adapted to a coastal fjord region (Collins et al., 2009). The model reproduces the mixing due to turbulent vertical velocities of unresolved eddies and numerically calculates horizontal velocities, temperature, salinity, and diffusivities in both space and time. To carefully determine the effects of surface mixing on phytoplankton growth in highly stratified regions, fine vertical resolution is needed. In Rivers Inlet, the initial spring bloom occurs near the surface with the chlorophyll maximum occurring above 3 m. This resolution is easily achieved in a one-dimensional model, but would be very difficult in a two- or three-dimensional model. Here, we use a vertical resolution of 0.25 m to resolve the near-surface mixing processes and to correctly depict the spring phytoplankton bloom. The necessary two-dimensional physical processes (estuarine circulation, baroclinic pressure gradients) are parameterized (Collins et al., 2009). Observations collected at Rivers Inlet are used to tune the one-dimensional model to correctly represent the Rivers Inlet physical processes.

In Rivers Inlet, the spring bloom is initiated as light limitation is lifted and is terminated by nitrate exhaustion. The first organisms to bloom are diatoms, primarily Thalassiosira spp. (David Cassis, University of British Columbia, pers. comm. and Shiller, 2012). Our model is thus limited to one nutrient type (nitrate) and one phytoplankton class (diatom). This model is fast and simple, allowing for both an interpolative mode for the years with observed data and an extrapolative mode for years where data do not exist, and for hypothetical forcings. This will allow us to be able to vary the physical forcings (windspeed and direction, cloud coverage, and freshwater flux) to determine the individual impacts on the timing of the spring bloom. The following sections will discuss the observational data needed to initialize and force the model, details of the physical model and biological model (Data and methods), model results (Results), discussion of the physical processes affecting the initiation of the spring phytoplankton bloom in Rivers Inlet, and the significance to other fjord regions (Discussion).

Data and methods

Data

The field data used in this research were collected for the Rivers Inlet Ecosystem Study during the early spring and summer of 2007, 2008, and 2009. Sampling took place approximately twice monthly between February and August (Hodal, 2011; Tommasi et al., 2013a). The station chosen to be modelled here is DF02 (Figure 1). In 2007 and 2008, daily sampling of temperature, salinity, photosynthetically active radiation (PAR), and chlorophyll fluorescence took place at the Florence Daily site located near DF02 down to 30 m using a Hydrolab DS5X sonde. During the 2008 season, the Hydrolab sonde broke down several times and an SBE 25 Sealogger CTD was used instead (Wolfe, 2010; Tommasi et al., 2013a). In 2009, the daily profiles were taken using a RBR XR-620 CTD. During the twice-monthly cruises, profiles of temperature, salinity, chlorophyll fluorescence, and PAR were taken using the SBE25 CTD. Nitrate, phytoplankton, and chlorophyll a were sampled using Niskin bottles at depths 0, 5, 10, and 30 m. The chlorophyll fluorescence measurements were calibrated to 0.7 µm filtered chlorophyll a.

Meteorological field data were obtained from two weather stations that were mounted near DF02 (Figure 1; R. Routledge, pers. comm., 2009; Hodal, 2011) and from the Port Hardy meteorological station (Environment Canada, 2009). Laska weather station replaced Ethel weather station in 2009 at the same location. Hourly values of windspeed and direction were collected from Ethel and Laska station during March and April 2008 and 2009. Data collected from the Port Hardy meteorological station included hourly values of humidity, cloud fraction, air temperature, and windspeed and direction. Daily Wannock River discharge data were obtained from the Water Survey of Canada (Environment Canada, 2006).

SOG model

The SOG model is a vertical, one-dimensional coupled biophysical model for the upper layers of a deep estuary. It consists of a mixing layer model (Large et al., 1994), parameterized estuarine processes (Collins et al., 2009), and a simple nutrient–phytoplankton biological model (Collins et al., 2009). It successfully predicts the timing of the spring phytoplankton bloom in the Strait of Georgia, British Columbia, Canada (Collins et al., 2009).

The model can be considered to consist of cylinder situated at the surface of the Inlet and extending down to 40 m depth, well above the total depth of the fjord (Figure 2). Within the cylinder, physical and biological quantities are assumed to uniform in the horizontal but vary in depth. Physically, the model solves for the non-estuarine flow considering time dependence, advection by the estuarine flow, the Coriolis force, pressure gradients, and mixing. It is forced by surface processes of heating (incoming solar radiation, outgoing longwave radiation, latent heat, heat conduction, and windstress). The estuarine flow is parameterized based on the river flow (more below) and it advects the temperature, salinity, velocities, and biological quantities. Near the surface, the flow entering the upstream side of the cylinder is more river-like (less saline, low in nutrients and chlorophyll) than the downstream side of the cylinder (Figure 2). This effect appears as dilution term in the salinity equation. In addition, the estuarine flow causes a vertical entrainment velocity which advects all properties. The vertical flow in the bottom of the cylinder requires boundary condition of all quantities at 40 m.

Figure 2.

Sketch showing the model as a cylinder at the surface in the Inlet. Solid near-horizontal line shows the surface. Dashed line shows the mixing layer depth; in the model domain, this is determined by the KPP model (Large et al., 1994). Not shown is the (i) mixing occurring throughout the cylinder, also determined by the KPP model and represented by $\mathcal{V}$ terms in Equation (1) and (ii) the Coriolis force, also determined by the KPP model. The physical model is forced at the surface using standard ocean model forcing parameterization (Collins et al., 2009). The solar heating penetrates the upper layers of the domain and appears in the second term on the LHS of Equation (1d). In the near-surface layers of the domain, the estuarine flux is from the river to the ocean and this flux is less saline going into the model domain than leaving. This appears in the model as a dilution term, the second term on the LHS of Equation (1c). The estuarine circulation also causes a vertical flow throughout the model domain w_e and this appears in all the equations (1). The bottom boundary conditions are the fluxes of heat and salt due to the vertical entrainment flux.

The physical model is tuned to provide realistic parameters to describe the physical dynamics of the Rivers Inlet system. The tuned physics model is coupled with the biology model which is then tuned independently to correctly match the timing of the spring phytoplankton bloom for the 3 years for which this timing is known.

Physical model

The physical model incorporates many additions to the original KPP mixing model to adapt it to a coastal region (Collins et al., 2009). Further modifications were made to use the model in the Rivers Inlet basin (Wolfe, 2010).

At this point, it is useful to distinguish between variables estimated from observations and tuned variables. There are many parameterizations used in the model in which a variable is determined from observational data of some kind and not tuned thereafter. A tuned variable, on the other hand, is one that is used in the model and is adjusted to produce the most favourable model output. Parameterizations of variables estimated from observations are discussed in the following subsections, while tuning variables are described in the Physical model tuning and Biological model tuning sections.

The physical model solves the following equations for across-fjord and along-fjord non-tidal, non-estuarine velocity, salinity, and temperature, respectively:

$$\frac{\mathrm{\partial }u}{\mathrm{\partial }t}-fv=\underset{\_}{-\frac{1}{{\rho }_o}\frac{\mathrm{\partial }p}{\mathrm{\partial }x}-\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{e}u)}+\mathcal{V}(u),$$

(1a)

$$\frac{\mathrm{\partial }v}{\mathrm{\partial }t}+fu=\underset{\_}{-\frac{1}{{\rho }_o}\frac{\mathrm{\partial }p}{\mathrm{\partial }y}-\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{e}v)}+\mathcal{V}(v),$$

(1b)

$$\frac{\mathrm{\partial }S}{\mathrm{\partial }t}=\underset{\_}{-\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{e}S)+J_s\exp (\gamma z)}+\mathcal{V}(S),$$

(1c)

$$\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\underset{\_}{-\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{e}T)}+\frac{I_{\mathrm{TOTAL}}(z)}{c_p\rho (z)}+\mathcal{V}(T),$$

(1d)

where u and v are the across- and along-fjord velocities, S is the salinity, T the temperature, t time, x and y are across- and along-fjord distance, z the vertical distance upward, starting at the sea surface (so all z values are negative), f the Coriolis parameter, $\mathcal{V}$ the mixing prescribed by the KPP model, including non-local mixing (Large et al., 1994), p the baroclinic pressure, ρ the density of the water, ρ_o a constant reference density, c_p the specific heat at constant pressure, w_e a vertical upwelling flux, J_s an added freshwater flux, 1/γ is the scale depth of the freshwater flux, and I_TOTAL the incoming shortwave radiation. Derivation of the momentum equations is given in Appendix 1. Surface boundary conditions for the momentum equations include the windstress and those for the temperature the longwave radiation balance and the sensible and latent heat flux.

The underlined terms are those additional to the standard KPP model. The horizontal pressure gradient terms account for the finite size the Inlet (Collins et al., 2009). The impact of the estuarine circulation in the Inlet is seen in two types of terms. In the salinity equation, the term J_s exp(γz) gives the dilution effect of the down-inlet flow of the upper fresh layer. In all equations, the term including w_e is the impact of the upward entrainment flux. These estuarine terms need to be parameterized and the derivations are given below.

Estuarine flow

Freshwater is added into the surface layer based on a parameterization of total freshwater flow into the inlet. The surface salinity at station DF02 is observed to be correlated with the Wannock River flux. An empirical fit between the two is:

$$S_{\mathrm{surface}}=S_D\frac{\exp (-Q_W/{\alpha }_1)+\beta \exp (-Q_W/{\alpha }_2)}{{\gamma }_1+\exp (-Q_W/{\alpha }_1)+\beta \exp (-Q_W/{\alpha }_2)},$$

(2)

where Q_W is the Wannock River flux, S_D = 31.8, the observed average salinity between 40 and 50 m in Rivers Inlet, α₁ = 80.0 m³ s⁻¹, α₂ = 1500.0 m³ s⁻¹, γ₁ = 0.04, β = 0.01. The empirical equation was chosen to agree with physical requirements such that at low river discharge values, the surface salinity approached the bottom boundary salinity value of 31.8 and with very high river discharge, the surface salinity goes to zero (Wolfe, 2010). The surface salinities of the model are then forced to this fit using the following expression for the dilution rate in the prognostic equation for salinity [Equation (1c)]:

$$J_s\exp (\gamma z)=F_w{Q}_W{S}_o{\left(\frac{Q_W}{\overline{Q}}\right)}^{{\sigma }_s}\exp \left(\frac{z}{a_h{h}_m}\right),$$

(3)

where F_w = 2.0 × 10⁻⁶ m⁻³, S_o the surface salinity at the previous time-step, a_h = 3.5, $\overline{Q}=7000$ m³ s⁻¹, σ_s = 1.38, and h_m the mixing layer depth determined by the KPP model. Again the shape of the function was chosen to replicate the underlying physical processes. The value of a_h was tuned to best match the model halocline depth to the observed halocline depth. The value of F_W was chosen to best match the mean model salinity to the mean observed salinity [Equation (2)] and the coefficient σ_s was chosen to make the relationship between modelled and observed salinity as linear as possible. The resulting correspondence between the modelled salinity and the derived fit [Equation (2)] has a Willmott index of agreement of 0.98 (Willmott et al., 1985) and root mean square error of 2.8 over a salinity range of 25.

Estuarine entrainment largely controls upwelling in a fjord and is also parameterized in the model. Following Collins et al. (2009), we use:

$$w_e=w_{\ast }\left[\frac{Q_W}{\overline{Q}}\exp \left(\frac{-Q_W}{\overline{Q}}\right)\right]\left[1-{\left(1-\frac{z}{2.5d}\right)}^2\right],$$

(4)

to relate the upwelling velocity to river discharge. Its structure was determined using a force balance between the estuarine pressure gradient and turbulent momentum mixing (Collins et al., 2009).

The magnitude of the vertical velocity w* = 1.08 × 10⁻⁴ m s⁻¹ as a function of Wannock River flow was determined using the Knudsen relations (Dyer, 1973) based on salinity observations at DFO2 assuming a 16 m deep surface layer. The value d = 6.4 m is the average depth at which the cumulative freshwater reaches 67% of its total value and was calculated using salinity profiles from DF02. For the flow rates observed for the Wannock River (<3000 m³ s⁻¹), Equation (4) is nearly linear in Q_W; we include the exponential term to give model stability for hypothetical very large river flows. This vertical velocity advects temperature, salt, momentum, nitrate, and phytoplankton in the model. The vertical velocity is strongest at depth compared with the surface causing a vertical convergence of water which results in a horizontal removal (advective loss) of these properties throughout the depth of the modelled water column.

Wind-driven vertical advection

The Strait of Georgia was modelled as a closed basin system (Collins et al., 2009). However, since advection out of the mouth of Rivers Inlet due to outflow winds was thought to be an important loss mechanism, it was necessary to include a parameterization for this process. A one-sided open basin model was developed to correctly depict the biophysical dynamics occurring within a fjord during outflow winds. During an outflow wind event, the mouth of the fjord is assumed to have an outflow velocity equal to the velocity seen at the central, modelled station. This causes the surface layer to flow out of the inlet, resulting in deeper water moving upwards to replace it. Upwelling occurs requiring the additional parameterization of a wind-driven vertical velocity which was added to the entrainment velocity, w_e. Vertical advection is thus being driven by both estuarine circulation and wind. Details of the one-sided open basin system can be found in Appendix 2.

Internal mixing

Rivers Inlet is much narrower than the Strait of Georgia. In numerical calculations, the narrowness of the inlet resulted in a baroclinic velocity structure with high vertical mode numbers, leading to step-like mixing which is not seen in the observations. To remove this unrealistic behaviour, the model was forced to have stronger implicit mixing by smoothing the diffusivity calculated by the model over 5 grid points rather than the 3 grid points used in Collins et al. (2009).

Light

The SOG physical model incorporates light in three parts: a cloud model (which determines the amount of light penetrating the atmosphere at a given time and hitting the water surface), an albedo (which controls how much light enters the water column), and the parameterization for the depth profiles of PAR, I_par, and the total light spectrum that is available for heating, I_total. The total non-reflected light entering the water column, I_total (0 m), is calculated by the geometry of the solar angle for a given time and latitude and the cloud filtering (Collins et al., 2009). The profiles of I_total with depth are based on the Jerlov (1976) classification (Collins et al., 2009).

To specify the light conditions in the water column, a photosynthetically active irradiance profile must be determined through a parameterization of the attenuation coefficient, K_par. Examination of I_par profiles taken at the daily station from February to June 2009 showed a positive correlation between the attenuation coefficient and the phytoplankton biomass in the water column (based on chlorophyll concentrations), but no significant correlation with river discharge. River turbidity does not reach down inlet to DFO2 until after the freshet begins, which occurs well after the spring bloom. Data from other years were unavailable. For better statistical robustness, we did an empirical fit of eK = exp(−K_pardz) as a function of depth with dz = 1 m rather than directly fitting K_par. The empirical relationship between the measured light and measured chlorophyll, P_c, was found to be

$$K_{\mathrm{par}}={\alpha }_l+{\beta }_l{\overline{P}}^{0.665}+{\theta }_l\exp \left(\frac{z}{d_l}\right),$$

(5)

where $\overline{P}=P_c/(1\text{mg}\text{Chl}{\text{m}}^{-3})$⁠, the exponent on $\overline{P}$ follows Ménesguen et al. (1995), α₁ = 0.1709 m⁻¹, β₁ = 0.02 m⁻¹, θ₁ = 2.53 m⁻¹, and d₁ = 0.53 m.

Wind

Windspeed and direction measurements come from Port Hardy Airport and the Ethel and Laska weather stations in Rivers Inlet (Environment Canada, 2009; R. Routledge, pers. comm., 2009; Hodal, 2011, respectively). Sources of historical meteorological data are not available for Rivers Inlet. In 2007, 2008, and 2009, winds were measured in Rivers Inlet. To patch gaps in the wind data, and for the wind data in 2006, an extrapolation based on Port Hardy wind was used. To represent the local wind direction in Rivers Inlet, a regression tree was developed using a comparison between the Port Hardy Airport and the Laska Meteorological Station wind directions (Wolfe, 2010).

Windspeed from Port Hardy is used as a proxy for the windspeed in Rivers Inlet. The linear relationship is:

$$V_{RI}=m{V}_P,$$

(6)

where V_RI is the windspeed at Laska weather station, V_P is the windspeed at Port Hardy meteorological station, and m = 0.877. Windspeed and direction are transformed into meridional and zonal components and rotated to align with the major axis of Rivers Inlet. The windstresses were calculated from these velocities following Large and Pond (1982).

Bottom boundaries

Temperatures and salinities from the 2008/2009 Rivers Inlet Ecosystem Study (RIES) cruises were used to determine an annual fit for the bottom boundary inputs. Measurements were averaged between 40 and 50 m and an empirical fit was determined using a constant and seasonal component. For temperature:

$$T_B=y_1+y_2\sin (\omega t+y_3),$$

(7)

where T_B is the temperature at 40 m, y₁ = 7.63°C, y₂ = 0.63°C, y₃ = 3.04, ω = 2π/365.25 d⁻¹, and t is the time in yearday. For salinity:

$$S_B=y_1+y_2\sin (\omega t+y_3),$$

(8)

where S_B is the salinity at 40 m, y₁ = 31.66, y₂ = 0.46, and y₃ = 4.51.

Nitrate concentrations at 40 m were considered constant at 21 µM consistent with observations seen at multiple stations throughout Rivers Inlet (Wolfe, 2010). The value was estimated from bottle samples collected during RIES, supplemented with data from CCGS J.P Tully cruises (M. Robert, unpublished data, 2005–2007. IOS/OSD Data Archive: DFO). Phytoplankton at 40 m was set to zero.

Biological model

The primary objective of the biological component of the SOG model, when coupled with the physical model, is to predict the timing of the spring phytoplankton bloom. A simple model is used to determine the growth rate of phytoplankton. All units in the biological model are nitrogen atom concentrations. Following Collins et al. (2009), the phytoplankton and nitrate equations are

$$\begin{array}{ll}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}& =\left\{min[u_c(I_{\mathrm{PAR}}),u_n(N)]-R_m-\mathrm{{\rm Y}}\frac{Z(t,z)}{{\kappa }_P+P}\right\}\\ & \times P+\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{s}P)-\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{e}P)+\mathcal{V}(P)\end{array},$$

(9a)

$$\frac{\mathrm{\partial }N}{\mathrm{\partial }t}=-\{min[u_c(I_{\mathrm{PAR}}),u_n(N)]\}P-\frac{\mathrm{\partial }}{\mathrm{\partial }z}(w_{e}N)+\mathcal{V}(N),$$

(9b)

where P is the phytoplankton concentration, N the nitrate concentration, u_c the growth based on light limitation, u_n the growth based on nitrate limitation, R_m the rate of natural mortality, Y the maximum ingestion of phytoplankton by zooplankton, Z the zooplankton concentration, K_P the half-saturation for zooplankton grazing, and w_s the sinking rate of the phytoplankton (positive downward). The last two terms in each equation are the estuarine advective loss and the turbulent diffusion determined by the physical model.

As light limitation is lifted, a spring bloom is initiated which is then terminated by nitrate exhaustion, the limiting nutrient in Rivers Inlet (Shiller, 2012). Only one size class of phytoplankton is used (P nominally micro-phytoplankton, >20 µm) and only nitrate N is modelled as a nutrient. During the spring bloom, over 99% of the total phytoplankton abundance is diatoms. Since Thalassiosira spp. is the first dominant component of the bloom, we model only this species. While ammonium is being used by phytoplankton, the amount of ammonium available at any time is too small to sustain the biomass of the spring bloom for longer than several hours, and its inclusion in the model was shown to make no substantial changes to the spring bloom timing elsewhere (Collins et al., 2009).

Phytoplankton concentrations are affected by growth rates, natural mortality, grazing, sinking, and physical processes, the left-hand side terms in Equation (9a), respectively. Growth rates are based on light and nutrient limitation. Following Collins et al. (2009), light-limited growth is given by

$$u_c=R_{\max }(T)\left[1-\exp \left(\frac{-I_{\mathrm{PAR}}}{0.67I_{\mathrm{opt}}}\right)\right]\left[1.8\exp \left(\frac{-I_{\mathrm{PAR}}}{2.7I_{\mathrm{opt}}}\right)\right],$$

(10)

where I_opt = 38.4 W m⁻² is the optimum light intensity and R_max(T) the maximum growth rate.

When light is not limiting, phytoplankton growth depends on nutrient availability

$$u_n=R_{\max }(T)\frac{N}{\kappa +N},$$

(11)

where K = 2 µM is the half-saturation constant for nitrate (N).

Growth rates, grazing, and mortality are dependent on temperature. Grazing is based on a tuned concentration of zooplankton that is held constant. Winter concentrations of zooplankton in Rivers Inlet are similar between years and do not change until after the bloom (Tommasi et al., 2013b). A density-dependent grazing rate was used to describe the functional response of zooplankton grazing at very low chlorophyll levels. The grazing threshold is defined as the concentration of prey, below which, the predator stops feeding (Leising et al., 2003; Cugier et al., 2005). A value of 0.09 mg Chl m⁻³ was used for the chlorophyll predation threshold which lies between values used by Wroblewski et al. (1988) and Cugier et al. (2005). Diatoms are negatively buoyant, particularly when nutrient-limited (Waite et al., 1992). Following Cugier et al. (2005), nutrient limitation is calculated as

$$f_N=\frac{N}{\kappa +N},$$

(12)

and then sinking speed is

$$w_s=W_1{f}_N^{0.2}+W_2(1-f_N^{0.2}),$$

(13)

where W₁ = 0.5 and W₂ = 1.2 m d⁻¹ are the sinking rates for nutrient replete and nutrient-depleted diatoms, respectively (Cugier et al., 2005). Sensitivity to these rates are weak, changing W₁ by 0.25 m d⁻¹ affects the timing of the spring bloom by <4 d and changing W₂ by the same amount has no effect.

The physical processes include loss due to entrainment and turbulent diffusion determined by the physical model. Advective terms due to horizontal variations in phytoplankton concentration are not included in the model. The timing of the spring bloom is nearly spatially uniform in the inlet with differences in bloom timing at different sites a maximum of 3 d apart (Tommasi, 2008). Table 1 gives values for all biological parameters modified from Collins et al. (2009).

Table 1.

Biological model parameters^a.

Symbol	Parameter	Value	Source
Rmax(Tref)	Maximum growth rate (24 h)	2.2 d−1	Tuned, <Hitchcock (1980), >Durbin (1974)
Tref	Reference temperature	10°
Prth	Predation threshold concentration	0.05	<Cugier et al. (2005), >Wroblewski et al. (1988)
RM(Tref)	Mortality	0.075 d−1	<Spitz et al. (2003), >Denman and Peña (2002)
	Maximum temperature for growth	18°	Durbin (1974)
	Range of dec. due to temp.	8°
Kp	Half-saturation for zooplankton grazing	0.2 µM N	<Denman and Peña (2002)
Y(Tref)	Maximum ingestion, mesozoo	0.6 d−1	Cugier et al. (2005)
Zw	Zooplankton concentration	0.089 µM N	Tuned
	Chl m3 to µM N conversion	1.7:1	M. Maldonado, per. comm., 2010

^aThe “<” and “>” symbols indicate whether the value is greater than (>) or less than (<) the referenced value.

Results

Physical model tuning

The physical model was first run independently of the biological model and compared with observations across all years for which data were available. The model was deemed useful when the modelled physical profiles resembled the observed physical profiles available for March–June (e.g. Figure 3). Surface salinity values, and depth and steepness of the halocline were quantitatively compared. To modify the physical model, internal mixing and freshwater parameters were tuned. The vertical advection velocity determined through w* was based on analytic calculations and was not tuned. The internal wave mixing values were set to those used in the Strait of Georgia (Collins et al., 2009). The dilution rate, through parameters F_w, magnitude, and σ_s, exponent in Equation (3), was tuned to produce surface salinities that matched the empirical relationship [Equation (2)]. The depth range in which freshwater was added [a_h in Equation (3)] was tuned to obtain a good match between the model and observed halocline depth.

Figure 3.

Example of a typical profile comparison between model and observations. In the left panel, for 12 April 2009, the shape of the temperature compares well. The right panel for the same day shows reasonable agreement for the salinity profiles. Density is determined primarily by salinity and so the strength and position of the halocline is dynamically important in the surface mixing. The average depth of the halocline over all profiles agrees with the data.

Biological model tuning

Once the physical model was tuned, its parameters were fixed. Next the biological model was tuned. Two biological parameters were tuned to provide the best match between the timing of the peak magnitude of the modelled spring bloom and the observed spring bloom for the years 2007–2009: maximum phytoplankton growth rate under optimal conditions (R_max) and the concentration of zooplankton (Z_w). The 0–3 m averaged phytoplankton concentrations were compared with daily chlorophyll fluorescence measurements for each model run and surface nitrate was compared with discrete bottle sampled nitrate at 0 m (Figure 4). Averaged 0–3 m phytoplankton concentration was used because the initial spring bloom is surface focused with the observed chlorophyll maximum occurring above 3 m in the water column. For the model, the spring bloom peak was defined as the date of the maximum 0–3 m averaged phytoplankton concentration within 4 d of the 0–3 m averaged nitrate concentrations going below 0.1 µM. A single value was found for each tuning parameter (Table 1) to produce optimal results for the timing of the spring bloom for the 3 years, 2007–2009 (Table 2).

Table 2.

Observed and modelled spring bloom dates.

Run ID	Start date	End date	Observed bloom	Modelled bloom
R6	1 September 2005	1 November 2006	29 March–12 April 2006	7 April 2006a
R7	8 October 2006	1 November 2007	28 April 2007	27 April 2007
R8	29 August 2007	1 November 2008	3 April 2008	3 April 2008
R9	22 September 2008	1 November 2009	18 April 2009	20 April 2009

^aThe modelled peak date for 2006 (*) uses modified wind directions as discussed in the Effect of winds section.

Figure 4.

Modelled 0–3 m averaged phytoplankton (solid line) and observed daily 0–3 m averaged chlorophyll (circles) for (a) 2007, (b) 2008, and (c) 2009. Surface modelled nitrate (solid line) and surface bottle sampled nitrate from biweekly cruises (x's connected by dashed line) for (d) 2007, (e) 2008, and (f) 2009. In the left panels, the observed spring bloom date is marked by a filled circle and the modelled spring bloom data are marked with an arrow.

Spring bloom timing

The mixing layer depth in Rivers Inlet responds to the winds and the river discharge (Figure 5). River flow in February–April was stronger in 2007 (226 m³ s⁻¹) compared with 2006 (95 m³ s⁻¹), 2008 (100 m³ s⁻¹), and 2009 (76 m³ s⁻¹) corresponding to a dampened mixing layer depth in 2007 (Figure 5). For each the year, the cube-root averaged cubed winds from January 1 to April 1 were 2006: 5.1 m s⁻¹, 2007: 4.6 m s⁻¹, 2008: 3.9 m s⁻¹, and 2009: 4.0 m s⁻¹. Mixing layer depth increases during strong winds. The mixing layer depth in 2009 was deepest (Figure 5) due to low river input and several outflow wind events. Strong outflow wind events (>7.4 m s⁻¹) were removed from the modelled 2006 winds (see the Effect of winds section). Summing the outflow wind cubed for strong outflow wind events and normalizing by (7.4 m s⁻¹)³ gives 48, 12, and 24 strong outflow hours for 2007, 2008, and 2009, respectively.

Figure 5.

Daily-averaged windspeed cubed (grey shaded), modelled daily-averaged mixing layer depth (shaded black), and Wannock River flow (solid line) for 1 February–1 November for each year. Stronger winds at similar river flow cause deeper mixed layers; compare April 2008 with April 2006. Stronger river flow enhances the stratification and dampens the effect of the winds; compare March 2007 with March 2009.

The interannual variation in the timing of the spring phytoplankton bloom is predicted successfully with the modelled peaks occurring within ±2 d of the observed blooms (Figure 4, Table 2). Considerable interannual variability was seen in the timing of the spring bloom at DF02. The earliest bloom occurs in 2008, with observed 0–3 m averaged chlorophyll values peaking on 3 April. The latest bloom occurs in 2007, on 28 April (Figure 4). Blooms in 2008 and 2009 are observed to occur 15 d apart, although windspeed, river flow, cloud coverage, and winter zooplankton abundance are similar between years.

After the spring bloom, the model shows continued high chlorophyll (Figure 4) which clearly disagrees with the observations. The model uses a constant winter-level zooplankton concentration, whereas in Rivers Inlet, the zooplankton dynamically respond to the spring bloom (Tommasi et al., 2013b). A responsive zooplankton model is beyond the scope of this manuscript and thus more than a few weeks after the spring bloom, the results here can be considered an estimate of what would happen if zooplankton did not increase after the spring bloom.

To quantify the effect of different physical forcings on the initiation of the spring bloom, a series of 48 tests was done (Table 3). To isolate the impact of the different forcings, in each test, one physical variable [wind (w), Wannock River flow (r), or cloud coverage (c)] was replaced with data from another period, while all other variables remained at their previous values. For example, “R6wR7” is the 2006 spring bloom, initialized in September 2005 and run until November 2006 using data from 2005 to 2006 except wind data, which were replaced with winds from 2006 to 2007. Each test scenario was compared with its control year, for example, the timing of the R6wR7 bloom would be compared with the 2006 bloom date using 2006 meteorological data. Results of the sensitivity tests are given in Table 4 and are described below.

Table 3.

List of test runs to determine the impact of cloud coverage, wind, and freshwater on the timing of the spring bloom.

Run ID	Control years	Forcing varied	Value of changed forcing
R*wR6	R7, R8, R9	Wind	R6 wind
R*wR7	R6, R8, R9	Wind	R7 wind
R*wR8	R6, R7, R9	Wind	R8 wind
R*wR9	R6, R7, R8	Wind	R9 wind
R*rR6	R7, R8, R9	Freshwater	R6 freshwater
R*rR7	R6, R8, R9	Freshwater	R7 freshwater
R*rR8	R6, R7, R9	Freshwater	R8 freshwater
R*rR9	R6, R7, R8	Freshwater	R9 freshwater
R*cR6	R7, R8, R9	Clouds	R6 clouds
R*cR7	R6, R8, R9	Clouds	R7 clouds
R*cR8	R6, R7, R9	Clouds	R8 clouds
R*cR9	R6, R7, R8	Clouds	R9 clouds

Table 4.

Change in days of the timing of the spring bloom in the sensitivity runs.

	R6	R7	R8	R9
R*wR6	–	−14	3	−13
R*wR7	7	–	10	−2
R*wR8	1	−19	–	−12
R*wR9	7	−8	10	–
R*rR6	–	−4	4	−8
R*rR7	2	–	6	−3
R*rR8	−10	−13	–	−7
R*rR9	−3	−5	4	–
R*cR6	–	−4	0	1
R*cR7	2	–	3	5
R*cR8	0	−4	–	1
R*cR9	0	−9	−1	–

Run naming convention given in Table 3.

Effect of winds

The timing of the spring bloom is highly sensitive to the windspeed and direction (Table 4). Wind events were particularly strong and frequent in 2006 and 2007 (Figure 5). However, the sensitivity study shows that 2007 and 2009 winds caused delayed blooms and the winds in 2006 and 2008 caused advanced blooms. Compared with 2007 and 2009, 2008 had weaker winds with few outflow wind events (Figure 6); changing to 2008 winds for 2007 and 2009, advanced the bloom by 19 and 12 d, respectively. Compared with 2007 and 2009, 2006 had similar strength winds but was modelled with fewer outflow wind events (Figure 6); changing to 2006 winds for 2007 and 2009, advanced the bloom by 12 and 14 d, respectively.

Figure 6.

Daily averaged windspeed cubed (grey shading), modelled 0–3 m averaged phytoplankton (solid black line), and daily-averaged along inlet wind component with inflow winds positive and outflow winds negative (red line) for all 4 years. For 2006, the five strong wind events (e1–e5, Table 5) have been set as inflow wind events. Winds were strong in 2007 with a number of strong outflow events (b) and winds were weakest in 2008 with few strong outflow events (c); the bloom is latest in 2007 (b) and earliest in 2008 (c).

Sensitivity to wind direction

For 2006, wind measurements in Rivers Inlet were not available, and Port Hardy winds were extrapolated. However, extrapolating wind direction from Port Hardy to simulate conditions in Rivers Inlet was not precise and resulted in possible error in wind direction even for strong wind events. Five strong outflow wind events (cubed root mean wind-cubed >7.4 m s⁻¹) in the months January–April were extrapolated from the Port Hardy winds (Table 5). These strong winds grouped naturally into five events, each of <6 d duration. The original 2006 run, with all five events as outflow winds, resulted in a spring bloom occurring on 27 April 2006 (Figure 7). To explore the sensitivity of the model to wind direction, a test case model was run where these wind events were reversed. With all five events as inflow winds, the spring bloom occurred on 7 April 2006 (Figure 7), in agreement with the observations (Table 2).

Table 5.

List of five strong wind events in winter/spring 2006.

Wind event label	Date of wind event	Average windspeed cubed
e1	2 January–8 January 2006	1569 m3 s−3 for 6 d
e2	31 January 2006	702 m3 s−3 for 1 d
e3	25 February, 28 February–1 March 2006	881 m3 s−3 for 3 d
e4	2 April 2006	611 m3 s−3 for 1 d
e5	15 April–19 April 2006	591 m3 s−3 for 4 d

Figure 7.

Daily-averaged windspeed cubed (grey shaded), modelled daily-averaged mixing layer depth (black shaded), modelled 0–3 m averaged phytoplankton (black line), and 3-d averaged along inlet wind component with inflow winds positive and outflow winds negative (red line) for (a) 2006 using winds derived from Port Hardy and (b) 2006 where large outflow winds are reversed in direction (cubed winds >400 m³ s⁻³). Large outflow wind events occurring on 24 March and 15 April 2006 in the Port Hardy derived winds clearly affect the mixing layer depth. This figure is available in black and white in print and in colour at ICES Journal of Marine Science online.

To further explore sensitivity, all the identified wind events but one were changed in direction to inflow to determine the effect of single outflow wind events on the timing of the spring bloom (Table 6). Then test runs were done with multiple outflow wind events (Table 6). Reversing all five wind events to inflow winds produced a bloom that occurred on 7 April 2006, a 20 d shift. Reversing outflow wind events that occurred in January resulted in a 3 d shift in the spring bloom timing. Reversing outflow wind events that occurred in March resulted in a shift of 5 d, whereas reversing events occurring in April shifted the bloom by up to 7 d. A test was done to determine the effect of removing an outflow wind event compared with reversing the winds to an inflow wind event (Table 6). Removing an outflow wind event occurring in mid-April caused the spring bloom to occur on 19 April 2006, 1 d earlier than the effect of reversing the wind event to an inflow event. That is, reversing the wind event advanced the bloom by 7 d, whereas completely removing it advanced the bloom by 8 d.

Table 6.

Sensitivity of spring bloom date to reversing/removing outflow wind events.

Wind event(s) reversed/removed	Modelled date of bloom
None (5 outflow events)	27 April 2006
e2	24 April 2006
e3	22 April 2006
e4	24 April 2006
e5	20 April 2006
e3, e5	16 April 2006
e1, e3, e5	13 April 2006
e2, e3, e5	11 April 2006
e5 removed	19 April 2006
All (five inflow events)	7 April 2006

Effect of Wannock river flow

The timing of the spring bloom is sensitive to freshwater flux. Average Wannock River discharge from February through April was stronger in 2007 (226 m³ s⁻¹) compared with 2006 (95 m³ s⁻¹), 2008 (100 m³ s⁻¹), and 2009 (76 m³ s⁻¹). This resulted in most modelled test runs computed with 2007 river data to experience a delayed bloom. However, the largest effect is with the 2008 river flow. The March 2008 river flow was moderate at 107 m³ s⁻¹, whereas the March 2007 river flow was strong at 233 m³ s⁻¹ and the March river flows in 2006 and 2009 were low at 65 and 71 m³ s⁻¹, respectively. Changing to 2008 river flow advanced all blooms (7–13 d). Switching to the strong 2007 river flow delayed the 2008 bloom (by 6 d) but only changed the bloom slightly when switched with the weak 2006 or 2009 river flows (2–3 d).

Wannock River flux does not respond directly to storms as it is moderated by the lake. Discharge varies over weeks not days. This slowly varying river flux has two opposite effects on the bloom timing. First, higher river discharge causes the water column to become more stratified, reducing the depth of the mixing layer for a given windforcing. For example, compare the depth of the mixing layer during the storm about 21 March 2006 during weak river flow, with its depth during the storm about 5 February 2007 during moderate river flow, with its depth during the storm about 10 June 2007 during high river flow (Figure 5). A shallower mixing layer depth means more light for the phytoplankton leading to an earlier bloom.

Second, large river discharge also causes higher estuarine entrainment directly in the model (4) and is shown by a box model inverting the observed salinity structure in Rivers Inlet (Hodal, 2011). Higher estuarine entrainment causes a larger advective loss term for phytoplankton by increasing the near-surface outflux (Figure 2). This loss delays the spring bloom.

A sensitivity test was run to determine which river effect was dominant in Rivers Inlet. The year 2006 was used as the run year for the test case. Constant river discharge was used between 1 January 2006 and 30 April 2006, while observed river discharge values were used the rest of the year. Sixteen test cases were run, starting with a constant discharge of 50 m³ s⁻¹ and increasing in increments of 10 m³ s⁻¹, to finish with a constant discharge of 200 m³ s⁻¹. A moderate river discharge of 70–140 m³ s⁻¹ resulted in a spring bloom that occurred on 7 April 2006 (year day 97). At low river discharge (50 m³ s⁻¹), the spring bloom was delayed by 3 d due to a deeper mixing layer depth. High river discharge values (200 m³ s⁻¹) resulted in a 4 d delay of the spring bloom due to high advection rates (Figure 8).

Figure 8.

A sensitivity test to determine which river effect was dominant in Rivers Inlet. (Top panel) The entrainment flux (at 40 m) caused by the river flow. Entrainment flux increasing the near-surface outflow causing advective loss of phytoplankton (Figure 2). (Middle panel) The average mixing layer depth between 1 January and 30 April 2006. (Lower panel) The resulting spring bloom date. Large river discharge (200 m³ s⁻¹) results in stronger entrainment leading to a larger advective loss term for phytoplankton, delaying the spring bloom. At low river discharge (50 m³ s⁻¹), the spring bloom is delayed due to a deeper mixing layer depth.

Effect of cloud coverage

Interannual differences in cloud coverage have a small effect on the timing of the spring bloom in Rivers Inlet. Cloud coverage is very consistent between years during the 2 months preceding the spring bloom (February and March). Average cloud coverage during February and March varies <10% between 2006 and 2010. The largest effect was seen by changing the clouds between 2007 and 2009. Using 2007 clouds resulted in a 5 d delay in the timing of the 2009 spring bloom, whereas using 2009 clouds advanced the timing of the 2007 spring bloom by 9 d. Cloud coverage from any of the other years resulted in an advance or delay of 4 d or less of the timing of the spring bloom.

Sensitivity to biological parameters

While only two biological parameters were used to tune the model [R_max (T_ref) and Z_w], three parameters were examined to determine the effects on the timing of the spring bloom: maximum growth rate [R_max (T_ref)], mortality [R_M(T_ref)], and concentration of zooplankton (Z_w). The timing of the spring bloom is highly sensitive to growth rate and mortality and less so to concentration of zooplankton. Increasing the maximum growth rate by 33% from 2.16 to 2.8 d⁻¹ results in a bloom that occurs 22 d earlier on average. Decreasing the mortality rate by 33% from 0.075 to 0.05 d⁻¹ results in a bloom that occurs on average 7 d earlier. Decreasing the concentration of zooplankton by 33% from 0.089 to 0.06 results in a bloom that occurs on average 3 d earlier.

The spread of the spring bloom is defined as the time difference of the spring bloom between different years. Growth rate and mortality of phytoplankton have a larger effect on phytoplankton concentrations during spring and summer than during winter. This results in changes in the spread of the spring bloom date when these parameters are varied. Increasing the maximum growth rate by 33% resulted in an average decrease in spread of 10 d, while decreasing the mortality rate by 33% resulted in an average decrease in spread of 11 d. The zooplankton concentration has a constant effect on the winter and spring concentrations of phytoplankton. This results in zooplankton concentration not having a significant effect on the spread between the spring bloom dates. Growth rate and mortality of phytoplankton have a similar effect on the timing of the spring bloom; thus, only one parameter is needed for tuning the model. The mortality was set to 0.075 d⁻¹ at 10°C and the growth rate was tuned as outlined in the Biological model tuning section.

Many of the biological parameters are not known to any precision and a different set could give a similar bloom timing. Considering this fact, Collins et al. (2009), using a similar model for the Strait of Georgia, investigated the bloom timing variation due to winds, river flow, and clouds, for three different sets of biological parameters. No qualitative differences were seen in the results and the same relative dependence on the forcing factors was robust to even large changes in the biological parameters.

Discussion

How do different physical forcings affect the timing of the spring bloom?

Windspeed and direction are the primary physical factors controlling the timing of the spring bloom. In winter, strong winds increase the mixing layer depth (Figure 5), limiting phytoplankton growth (Figure 6). Increasing the mixing layer depth causes the phytoplankton to be mixed below their critical depth at which there is enough light to grow. Many large wind events are seen to decrease the depth-averaged chlorophyll concentrations, and windspeed and direction in the 2 months before the spring bloom have a significant impact on the timing of the spring bloom.

However, our results also show that the direction of large wind events can amplify the mixing layer depth. If the wind direction is directed down inlet towards the mouth, the freshwater layer is advected out of the fjord which is conducive to deeper mixing. Large wind events that are directed up inlet towards the head do not have this effect as the freshwater layer is not lost from the inlet. For example, the late February, early March outflow wind event in 2006, led to removal of much of the surface layer which then allowed deep mixing by the early March storm (Figure 7a). Whereas with the late February, early March winds reversed, the surface layer was not removed and mixing due to the early March storm is much reduced (Figure 7b). However, outflow wind events will also result in flushing events that advect phytoplankton populations away from the basin further delaying the timing of the bloom. A single outflow wind event within 2 weeks of the spring bloom causes a 7 d delay in the bloom timing (Table 6).

Freshwater flux is also an important control. Higher river flux increases the upwelling advection while shallowing the mixing layer depth. Low river discharge causes the opposite effect: less advection but a deeper mixing layer depth. Both effects were seen among the 4 years modelled; 2006 and 2009 river flows in March were low, allowing deep mixing layers, whereas in 2007, river flow in March was high, leading to strong advection. River flows in March 2008 were moderate and switching to 2008 river flows advanced all spring blooms (Table 4).

Even given that light limitation is an important process, interannual cloud coverage has a smaller effect on the timing of the spring bloom in Rivers Inlet. Cloud coverage is very consistent between years in Rivers Inlet with an average of 80% cloud coverage in February and March. Light is considered an important controller of phytoplankton biomass accumulation (Reynolds, 2006) and is seen as a controlling factor in many locations including the Strait of Georgia (Collins et al., 2009) and northerly fjords in Norway (Eilertsen et al., 1995). If climatic changes were to effect the average winter and early spring cloud distribution over Rivers Inlet, cloud cover would have a much larger potential influence on the timing of the spring bloom.

When acting through the mixing layer depth, the impacts of wind, wind direction, and freshwater flux on spring bloom timing are consistent with the Sverdrup (1953) hypothesis. However, the advective impacts of wind-driven flow and estuarine flow mean that phytoplankton loss are also important in the interannual variation of spring bloom timing in this system. Advective loss has also been shown to be a driver of spring bloom timing in, for example, Otsuchi Bay, Japan (Furuya et al., 1993), and the San Francisco Estuary, USA (Lucas et al., 1998). In the former case, the advective loss was horizontal, as it is in Rivers Inlet. In the latter case, the advective loss was vertical due to diatom sinking.

Wind-driven advection

Blooms in 2008 and 2009 were observed to occur 15 d apart. The mean cubed root averaged windspeed was similar (71 and 68 m⁻³s⁻³) through the 1 January–14 March period for the 2 years. However, the mean cubed root average wind vector (taking into account the direction of the wind) for the same period is quite different (45 and 30 m⁻³s⁻³, respectively); 2009 had more outflow wind than 2008. In particular, it had a strong outflow event in early March (Figure 6) which decreased the phytoplankton biomass as much as the larger wind events bracketing it (28 February and 16 March). Not including wind-driven advection out of the fjord led to bloom timings within a couple of days of each other. Thus, the difference in timing is due, in part, to wind-driven advection by outflow winds. Because of the high stratification seen in Rivers Inlet due to continual river input and the non-constricted opening at the mouth of the fjord, a large wind event directed down inlet causes rapid horizontal advection, flushing the phytoplankton population out of the inlet. The intensity of the flushing events depends on outflow wind timing and strength, strength of water column stratification, and how constricted the mouth opening is. Wind events large enough to reverse surface-layer flow were found to be a common occurrence in Knight Inlet (Baker and Pond, 1995). During 2–3 d strong wind events in Knight Inlet, speeds at 2 m depth were as large as 30 cm s⁻¹ and up to 5 cm s⁻¹ as deep as 20–25 m. Incorporating an open basin system in our Rivers Inlet model resulted in an additional loss term of phytoplankton due to the advection from wind which produced significant effects on the timing of the spring bloom and explained the difference between the bloom timing in 2008 and 2009.

An important result is the finding that single outflow wind events have significant effects on the timing of the spring bloom. A single outflow wind event can delay the timing of the spring bloom by 7 d. Removing an outflow wind event resulted in an 8 d earlier spring bloom, while the removal of an inflow wind event affected the spring bloom by less than a day. An outflow wind event will result in two mechanisms that delay the bloom timing: basin flushing and deeper mixing. An inflow wind event will not instigate basin flushing, but will result in deeper mixing. We would expect that removing a wind event of either direction would cause an earlier spring bloom due to a shallower mixing layer depth. It appears that in Rivers Inlet, while an inflow wind event does create a deeper mixing layer, the dominant effect on phytoplankton growth is determined by flow direction.

Significance

In the narrow valleys and fjords in British Columbia, significant winds will only occur as inflow or outflow (up or down valley; Jackson and Steyn, 1994). Outflow is the dominant direction in winter (Jackson and Steyn, 1992). Strong outflow winds, also known as gap or Squamish winds, occur when cold, dry air at low levels moves from the interior to the coast, dissecting the coastal mountains (Jackson, 1993). Outflow winds occur when a strong Arctic anticyclone develops a high pressure cold and stable air mass over the interior, east of the Coast Mountains. Differences in temperature and humidity between this air mass and the warm, moist coastal air create a large pressure gradient that is oriented perpendicular to the coast. This will cause strong low-level winds that develop in fjords and valleys. Strong outflow can also occur as a result of a deep cyclonic centre that approaches the coast. This will create a similar pressure gradient, driving outflow winds through the fjords. These deep cyclonic centres produce frequent rainfall and southwesterly winds during winter in the Pacific Northwest (Logerwell et al., 2003).

Inflow is the dominant direction in summer. The switch between outflow and inflow is correlated with the larger spatial scale change between winter southwesterly winds along the British Columbia coast and summer northeasterly winds along the coast. This change, termed the spring transition, can occur any time between March and June (Logerwell et al., 2003; Thomson and Hourston, 2009). In Rivers Inlet, large outflow wind events are not seen during summer (Figure 6). The closest long-term weather station to the north (Cathedral Point) also shows very large outflow wind events occurring only in winter and spring months. During summer (June–August), smaller inflow wind events occur at Cathedral Point. The timing of strong outflow wind events with windspeed >14 m s⁻¹ were calculated from 1995 to 2008 using data from Cathedral Point (Figure 9). The latest recorded outflow wind event occurred on 20 April in 1997. However, while Rivers Inlet and Cathedral Point both experience an increase in large outflow wind events caused by large-scale low pressure systems during winter, the timing of such events is not correlated. Therefore, while the spring transition date will likely affect the occurrence of outflow winds for all inlets along the coast of British Columbia, the stochastic nature of these wind events does not allow a non-local generalization of their timing. Unfortunately, this means it is difficult to hindcast bloom dates for Rivers Inlet in the same manner as was done in the Strait of Georgia (Allen and Wolfe, 2013).

Figure 9.

Outflow wind events occurring at Cathedral Point with windspeed values >14 m s⁻¹ (asterisk) and the spring transition date calculated by Thomson and Hourston (2009) (O and dashed line).

However, some generalization is still possible. The winds before the spring transition favour outflow within fjords which causes a delay in the timing of the spring bloom. Logerwell et al. (2003) indexed the spring transition date for a region off Oregon, USA, based on the first day where the value of the 10-d running average for upwelling is positive and the 10-d running average for sea level anomaly is negative. Based on this index, the mean date of the spring transition occurs on 6 April but can range from early February to early July (Figure 10). Thomson and Hourston (2009) determine the spring transition timing off of the south coast of British Columbia. The spring transitions showed strong similarities between the Oregon and British Colombia coastlines (Figure 10). A comparison between large outflow wind events occurring at Cathedral Point and the spring transition timing calculated by Thomson and Hourston (2009) shows that most strong outflow wind events occur before the spring transition (Figure 9). Large outflow wind events that occur later in spring are correlated with late spring transitions. Thus, late transition dates are more conducive to outflow winds occurring later in spring which will delay the onset of the bloom.

Figure 10.

Anomalies in the date of the physical spring transition from 1969 to 2008 for the OPI area (shaded bars) and spring transition dates for the west coast of Vancouver Island (Thomson and Hourston, 2009). Anomalies are based on an average date of 6 April (Logerwell et al., 2003).

It has been suggested that the timing of the spring bloom may be critical to the growth and survival of higher trophic levels (Cushing, 1974, 1975, 1982; Platt et al., 2003; Richardson, 2004; Ware and Thomson, 2005; Borstad et al., 2011). The sockeye salmon crash in Rivers Inlet occurred in the early 1990s, roughly coinciding with a period of late transition dates seen off Vancouver Island (Figure 10). It has been hypothesized that bottom-up effects control juvenile sockeye salmon growth in Rivers Inlet, and further, that high growth in the early marine phase is critical to at sea survival (McKinnell et al., 2001; Ainsworth et al., 2011). Changes in the Rivers Inlet spring phytoplankton bloom cause interannual variability in the timing, magnitude, and composition of the zooplankton bloom (Tommasi et al., 2013a, b). Match/mismatch of the phytoplankton bloom and the inlet zooplankton community may have resulted in less prey for juvenile salmon, thus diminishing marine survival and is a possible scenario contributing to the decline of the sockeye salmon in Rivers Inlet.

Acknowledgements

Rick Routledge initialized the Rivers Inlet project and led the collection of meteorological and oceanographic data in 2006 and 2007. The authors thank Evgeny Pakhomov for his project leadership from 2008. The excellent dataset used in this study would not have been collected without the efforts of many people including our sea-techs: Lora Pakhomova and Chris Payne, UBC graduate students: Asha Ajmani, Jade Shiller, Mike Hodal, and Desiree Tommasi, the Rivers Inlet family: Bachen Family, and the Western Bounty Crew. Special thanks go to Doug Latornell for technical support regarding the SOG model, Maite Maldonado for insight into the biological workings of Rivers Inlet, Jeremy Sklad for data analysis regarding the outflow wind events at Cathedral Point, and an anonymous reviewer for many constructive comments. Financial support for the Rivers Inlet project was provided by the Tula Foundation.

Appendix 1

Derivation of the physical equations

Starting with the rotating Navier–Stokes equations, where we have made the Boussinesq approximation and separated the flow into a mean and turbulent part, we have the following horizontal momentum equations [e.g. Pedlosky (1979)]:

$$\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }u}{\mathrm{\partial }y}+w\frac{\mathrm{\partial }u}{\mathrm{\partial }z}-fv=\frac{-1}{{\rho }_o}\frac{\mathrm{\partial }p}{\mathrm{\partial }x}+{\mathcal{V}}_u(u,v),$$

(14a)

$$\frac{\mathrm{\partial }v}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }v}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }v}{\mathrm{\partial }y}+w\frac{\mathrm{\partial }v}{\mathrm{\partial }z}+fu=\frac{-1}{{\rho }_o}\frac{\mathrm{\partial }p}{\mathrm{\partial }y}+{\mathcal{V}}_v(u,v),$$

(14b)

where (u,v,w) are the velocity components in the across-fjord, along-fjord, and vertical directions, f is the Coriolis parameter, ρ_o a constant reference density, p the pressure, and $\mathcal{V}$ the parameterization of the Reynolds' stresses as a function of the velocity. Furthermore, we assume the fluid is hydrostatic and incompressible [e.g. Pedlosky (1979)] so that:

$$\frac{\mathrm{\partial }p}{\mathrm{\partial }z}=-\rho g,$$

(15a)

$$\frac{\mathrm{\partial }u}{\mathrm{\partial }x}+\frac{\mathrm{\partial }v}{\mathrm{\partial }y}+\frac{\mathrm{\partial }w}{\mathrm{\partial }z}=0.$$

(15b)

In addition, we have three equations that describe (i) the conservation of salt, (ii) the heat and temperature, and (iii) the density as a function of salt and temperature [e.g. Pond and Pickard (1983)]. For the latter, we use the International Equation of State 1980 formulation (Pond and Pickard, 1983; their Appendix 3)

$$\frac{\mathrm{\partial }S}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }S}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }S}{\mathrm{\partial }y}+w\frac{\mathrm{\partial }S}{\mathrm{\partial }z}={\mathcal{V}}_S(u,v,S),$$

(16a)

$$\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }T}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }T}{\mathrm{\partial }y}+w\frac{\mathrm{\partial }T}{\mathrm{\partial }z}=\frac{dI}{dz}+{\mathcal{V}}_T(u,v,T),$$

(16b)

where S is salinity, T the temperature, I the incoming solar radiation, and $\mathcal{V}$ are the Reynolds' fluxes.

We can formally separate the circulation into three components, that due to the tide (v_T), that due to the estuarine circulation (v_e), and the remainder which is mainly driven by the wind (v_w). Momentum equations can be derived separately for each, but because the equations are not linear, the separation is not complete. For example, for the along-fjord direction, we can define:

$$v=v_T+v_e+v_w,$$

(17a)

$$\frac{\mathrm{\partial }{v}_T}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_T{v}_T)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{v}_T^2}{\mathrm{\partial }y}+\frac{\mathrm{\partial }(v_T{w}_T)}{\mathrm{\partial }z}+f{u}_T=\frac{-1}{{\rho }_o}\frac{\mathrm{\partial }{p}_T}{\mathrm{\partial }y}+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(\mathcal{K}(u,v)\frac{\mathrm{\partial }{v}_T}{\mathrm{\partial }z}\right),$$

(17b)

$$\frac{\mathrm{\partial }{v}_e}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_e{v}_e)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{v}_e^2}{\mathrm{\partial }y}+\frac{\mathrm{\partial }(v_e{w}_e)}{\mathrm{\partial }z}+f{u}_e=\frac{-1}{{\rho }_o}\frac{\mathrm{\partial }{p}_e}{\mathrm{\partial }y}+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(\mathcal{K}(u,v)\frac{\mathrm{\partial }{v}_e}{\mathrm{\partial }z}\right),$$

(17c)

$$\begin{array}{ll}\frac{\mathrm{\partial }{v}_W}{\mathrm{\partial }t}& +\frac{\mathrm{\partial }[u_W(v_T+v_e)]}{\mathrm{\partial }x}+\frac{\mathrm{\partial }(u{v}_W)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }(u_T{v}_e)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }(u_e{v}_T)}{\mathrm{\partial }x}\\ & +\frac{\mathrm{\partial }[v_W(v_T+v_e)]}{\mathrm{\partial }y}+\frac{\mathrm{\partial }(v{v}_W)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }(v_T{v}_e)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }(v_e{v}_T)}{\mathrm{\partial }y}\\ & +\frac{\mathrm{\partial }[w_W(v_T+v_e)]}{\mathrm{\partial }z}+\frac{\mathrm{\partial }(v_{W}w)}{\mathrm{\partial }z}+\frac{\mathrm{\partial }(v_e{w}_T)}{\mathrm{\partial }z}+\frac{\mathrm{\partial }(v_T{w}_e)}{\mathrm{\partial }z}\\ & +f{u}_W=\frac{-1}{{\rho }_o}\frac{\mathrm{\partial }{p}_W}{\mathrm{\partial }y}+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(\mathcal{K}(u,v)\frac{\mathrm{\partial }{v}_W}{\mathrm{\partial }z}\right),\end{array}$$

(17d)

where K is the effective eddy viscosity and conservation of volume [Equation (15b)] has been used to write the advection terms in flux form. Direct measurements of velocity in Rivers Inlet have not been made, but we can estimate the magnitude of these various velocities.

Tidal amplitudes are of order 3 m (Fisheries and Oceans, Canada, 2014). Thus, 30 km from the end of a fjord, in 300 m of water, tidal velocities would be on order of v_T ≈ 1.5 cm s⁻¹. For reference, this velocity gives a tidal excursion of 300 m. Tidal velocities are primarily back-and-forth along the fjord and the vertical velocity associated with them w_T ≈ 0.1 mm s⁻¹ simply raises and lowers the whole water column.

The estuarine circulation was calculated using a box model and the observed density structure (Hodal, 2011). During spring, before the freshet, outward flux in the upper layer (top 10 m) was found to be 700–1200 m³ s⁻¹. Given an approximate 3 km width to the fjord, horizontal estuarine velocities are on order of v_e ≈ 3 cm s⁻¹. Unlike tidal velocities, these small estuarine velocities have important consequences to the stratification. The horizontal flow is uni-directional bringing fresher water down the fjord. A vertical flow is induced which brings deeper water up into the upper layer (Hodal, 2011) w_e ≈ 0.2 µm s⁻¹.

Wind-induced surface flows in a nearby fjord have been measured as up to 30 cm s⁻¹ and rarely <5 cm s⁻¹ (Baker and Pond, 1995). Thus, v_W ≈ 10 cm s⁻¹ and by geometry v_W ≈ 10 µm s⁻¹.

Mixing approximately increases as velocity cubed [e.g. Pond and Pickard (1983)]. Thus, the ratio of mixing due to tides to mixing due to wind is 0.003 and between mixing due to estuarine circulation to mixing due to wind is 0.03. The combined effect of estuarine flow on wind mixing could increase or decrease mixing by 0.3. In the modelled, centre of the fjord, well away from the incoming river, mixing will be dominated by wind mixing. This is in agreement with many studies on estuarine systems [e.g. MacDonald and Geyer, 2004;,Hetland, 2005]

Thus, to calculate the mixing and Reynolds' stress terms, we need only calculate the wind-driven/remainder velocities (v_W). This brings our focus to Equation (17d) which is the tendency equation for the along-fjord velocity. In addition to the mixing term, there are a large number of advection terms. Most of these terms oscillate in time due to wind and/or tides. To determine how important these terms are we need to consider spatial gradients and time-scales. Winds in fjords vary rapidly due to variations in heating with the diurnal solar cycle. So a typical time-scale is 6 h. In the central part of the fjord, spatial scales are order 10 km. So the ratio of horizontal advection (using the dominate wind-driven velocities v_W) to time dependence is ∼0.2 and we neglect these terms in the momentum advection equation. Furthermore, the vertical gradient of the tidal velocities is negligible and the vertical tidal velocity simply raises and lowers the water column. As the horizontal estuarine circulation is small compared with the wind-driven flow, we also ignore its vertical advection. Thus, leaving us with:

$$\frac{\mathrm{\partial }{v}_W}{\mathrm{\partial }t}+\frac{\mathrm{\partial }[v_w(w_e+w_w)]}{\mathrm{\partial }z}+f{u}_W=\frac{-1}{\mathrm{\partial }{\rho }_o}\frac{\mathrm{\partial }{p}_W}{\mathrm{\partial }y}+\mathcal{V}(u_w,v_w),$$

(18)

which is equivalent to Equation (1b). Note that, if we assume the cross-fjord wind-driven velocities are the same order as the along-fjord velocities, which is definitely an overestimate, the Coriolis term is larger than the time-dependent term by a ratio of 2. It is actually much smaller as the cross-fjord velocities are small, but as this term is local and easy to handle in a one-dimensional model, we retain it. In addition, the vertical advection of the horizontal velocity is small (vertical lengthscale for upper layer ∼10 m), but for symmetry with the tracer equation, we retain this term as well. The pressure gradients due to the wind-driven flow are estimated from the density gradients generated by time-integrating the flow (Collins et al., 2009).

Appendix 2

Wind-driven vertical advection

For a semi-enclosed basin, such as the Strait of Georgia, the SOG code considers the basin to be a closed ellipse. The amount of water at each level moving in each direction is integrated and the thickening/thinning of the isopycnal layers is calculated (Collins et al., 2009). Baroclinic pressure gradients are calculated and this allows rotationally modified internal seiching to occur in basin.

In a fjord, such as Rivers Inlet, flows are predominately along axis. Flow towards the head of the inlet acts much like the flow in the closed basin described above. However, surface flow out of the inlet will leave the inlet and not cause back pressure gradients. This upper layer is lost to the inlet and deeper water moves upward. The background estuarine circulation is calculated separately. The strong surface flows considered here are primarily driven by the wind. Outward flows driven by the wind were observed to penetrate to ∼15 m depth in a nearby fjord (Baker and Pond, 1995). We use this depth to separate flow that easily leaves the fjord from deeper return-type flows.

For each depth-point in the model, the height of that isopycnal at each of the four directions east, west, north, and south which run across and along the axis of the eclipse is tracked. The density field at the east, west, north, and south extremes of the basin are integrated to calculate the baroclinic pressure field. At each step, the vertical mean is removed from the newly calculated velocities from the KPP mixing model to give the baroclinic currents. These velocities are added to time-integrated velocities along and across the fjord. For all directions, except the open end, the updated isopycnal distortions are calculated simply based on the time-integrated velocities (Collins et al., 2009).

For the open end, the isopycnals, above the depth of 15 m, are not allowed to tip downward. If the flow is outward in the surface layer, the isopycnals towards the mouth will be limited, whereas those towards the head will not be. Thus, these isopycnals will not be straight but kinked (Figure A1). If the velocity above the kinked isopycnal is outward towards the mouth of the inlet, a vertical velocity at the centre point is calculated as:

$$w_w={\int }_0^{\mathrm{bot}}\frac{v\kappa }{Ly/2}(-\text{d}z),$$

(19)

where bot = 40 m is the bottom of the model and v the outward velocity if the velocity is outward and zero otherwise. As a last step, the integrated velocity in each layer is adjusted to conserve mass given the divergence in the vertical flow.

Figure A1.

Sketch showing the upper part of the model water column during outward surface flow. The black circles show the calculated isopycnal positions at the head and mouth of the fjord. Solid lines connect them. The dotted line shows the extension of the isopycnal tilt and the definition of the kink, K.

Author notes

Handling editor: Rubao Ji