Kokanee

Data Preparation

The fish and downstream temperature data were provided by the Ktunaxa Nation. The discharge and elevation data were queried from the Columbia Basin Hydrological Database.

The data were cleaned and tidied using R version 4.0.5 (R Core Team 2020).

Length Cutoffs

Individuals were classified as fry (age-0), juvenile (age-1 and older subadults) or adult (sexually mature) based on the length cut-offs in Table 1.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The estimates were produced using JAGS (Plummer 2003). For additional information on Bayesian estimation the reader is referred to McElreath (2016).

Unless stated otherwise, the Bayesian analyses used weakly informative normal and half-normal prior distributions (Gelman, Simpson, and Betancourt 2017). The posterior distributions were estimated from 1500 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of 3 chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that the potential scale reduction factor \(\hat{R} \leq 1.05\) (Kery and Schaub 2011, 40) and the effective sample size (Brooks et al. 2011) \(\textrm{ESS} \geq 150\) for each of the monitored parameters (Kery and Schaub 2011, 61).

The parameters are summarised in terms of the point estimate, lower and upper 95% credible limits (CLs) and the surprisal s-value (Greenland 2019). The estimate is the median (50th percentile) of the MCMC samples while the 95% CLs are the 2.5th and 97.5th percentiles. The s-value can be considered a test of directionality. More specifically it indicates how surprising (in bits) it would be to discover that the true value of the parameter is in the opposite direction to the estimate. An s-value of 4.3 bits, which is equivalent to a p-value (Kery and Schaub 2011; Greenland and Poole 2013) of 0.05, indicates that the surprise would be equivalent to throwing 4.3 heads in a row. The requirement that non-essential explanatory variables have s-values \(\geq\) 4.3 bits provides a useful model selection heuristic (Kery and Schaub 2011).

The results are displayed graphically by plotting the modeled relationships between particular variables and the response(s) with the remaining variables held constant. In general, continuous and discrete fixed variables are held constant at their mean and first level values, respectively, while random variables are held constant at their typical values (expected values of the underlying hyperdistributions) (Kery and Schaub 2011, 77–82). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent or n-fold change in the response variable) with 95% credible intervals (CIs, Bradford, Korman, and Higgins 2005).

The analyses were implemented using R version 4.0.5 (R Core Team 2020) and the mbr family of packages.

Model Descriptions

Body Condition

The annual variation in condition (body weight when accounting for body length) was estimated from the boat and backpack electrofishing captures using a mass-length model (He et al. 2008).

Key assumptions of the condition model include:

• Weight varies with body length as an allometric relationship, i.e., \(W = \alpha L^{\beta}\).
• \(\alpha\) varies randomly with year.
• \(\alpha\) varies by period (pre versus post Mica 5 and 6).
• The residual variation in weight is log-normally distributed.

Preliminary analyses indicated that site and day of the year were not informative predictors of condition.

Relative Abundance

The annual variation in relative abundance was estimated from the boat count and catch data using an over-dispersed Poisson model (Kery and Schaub 2011). Lineal densities are by kilometre of river (as opposed to kilometre of bank).

Key assumptions of the relative abundance model include:

• Lineal count density varies by period.
• Lineal count density varies randomly with year.
• Lineal catch efficiency is a fixed multiplier of lineal count efficiency.
• Expected counts (and catches) are the product of the count (catch) density and the length of river (half the length of bank) sampled.
• Observed counts (and catches) are described by a Poisson-gamma distribution.

Preliminary analyses indicated that site and discharge were not informative predictors of the lineal count (or catch) density.

The model does not distinguish between the abundance and observer efficiency, i.e., it estimates the count which is the product of the two. As such it is necessary to assume that changes in observer efficiency by year are negligible in order to interpret the estimates as relative abundance.

Water Temperature

Model Templates

.model{
  bAlpha ~ dnorm(5, 2^-2)
  bBeta ~ dnorm(3, 2^-2)

  bAlphaPeriod[1] <- 0
  for(i in 2:nPeriod) {
    bAlphaPeriod[i] ~ dnorm(0, 2^-2)
  }

  sAlphaAnnual ~ dnorm(0, 2^-2) T(0,)
  for(i in 1:nAnnual) {
    bAlphaAnnual[i] ~ dnorm(0, sAlphaAnnual^-2)
  }

  sWeight ~ dnorm(0, 2^-2) T(0,)
  for (i in 1:length(Weight)) {
    eAlpha[i] <- bAlpha + bAlphaPeriod[Period[i]] + bAlphaAnnual[Annual[i]]
    eBeta[i] <- bBeta
    log(eWeight[i]) <- eAlpha[i] + eBeta[i] * Length[i]
    Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
  }

.model{
  bEfficiencyVisitType[1] <- 0
  for (i in 2:nVisitType) {
    bEfficiencyVisitType[i] ~ dnorm(0, 1^-2)
  }

  bDensity ~ dnorm(0, 5^-2)

  bDensityPeriod[1] <- 0
  for(i in 2:nPeriod) {
    bDensityPeriod[i] ~ dnorm(0, 2^-2)
  }

  sDensityAnnual ~ dnorm(0, 2^-2) T(0, )
  for(i in 1:nAnnual) {
    bDensityAnnual[i] ~ dnorm(0, sDensityAnnual^-2)
  }

  sDispersion ~ dnorm(0, 2^-2) T(0, )
  for (i in 1:length(Annual)) {
    log(eEfficiency[i]) <- bEfficiencyVisitType[VisitType[i]]
    log(eDensity[i]) <- bDensity + bDensityPeriod[Period[i]] + bDensityAnnual[Annual[i]]
    eAbundance[i] <- eDensity[i] * SiteLength[i] / 2
    eDispersion[i] ~ dgamma(1 / sDispersion^2, 1 / sDispersion^2)
    Count[i] ~ dpois(eAbundance[i] * eEfficiency[i] * eDispersion[i])
  }

.#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {

DATA_VECTOR(Tailrace);
DATA_VECTOR(Intake);
DATA_VECTOR(Forebay);
DATA_VECTOR(Air);
DATA_VECTOR(Unit1);
DATA_VECTOR(Unit2);
DATA_VECTOR(Unit3);
DATA_VECTOR(Unit4);
DATA_VECTOR(Unit5);
DATA_VECTOR(Unit6);
DATA_VECTOR(Spill);
DATA_IVECTOR(Consecutive);

PARAMETER(sDepthUnit);
PARAMETER_VECTOR(bDepthUnit);
PARAMETER(bDepthSpill0);
PARAMETER_VECTOR(bDepthTurbineUnit);
PARAMETER(bRelease);
PARAMETER(bReleaseAir);
PARAMETER(bRho);
PARAMETER(sTailrace);

vector<Type> eTurbineUnit1 = Tailrace;
vector<Type> eTurbineUnit2 = Tailrace;
vector<Type> eTurbineUnit3 = Tailrace;
vector<Type> eTurbineUnit4 = Tailrace;
vector<Type> eTurbineUnit5 = Tailrace;
vector<Type> eTurbineUnit6 = Tailrace;
vector<Type> eDischarge = Tailrace;
vector<Type> eDepthTurbineUnit1 = Tailrace;
vector<Type> eDepthTurbineUnit2 = Tailrace;
vector<Type> eDepthTurbineUnit3 = Tailrace;
vector<Type> eDepthTurbineUnit4 = Tailrace;
vector<Type> eDepthTurbineUnit5 = Tailrace;
vector<Type> eDepthTurbineUnit6 = Tailrace;
vector<Type> eDepthSpill = Tailrace;
vector<Type> eStratification = Tailrace;
vector<Type> eEffectTurbineUnit1 = Tailrace;
vector<Type> eEffectTurbineUnit2 = Tailrace;
vector<Type> eEffectTurbineUnit3 = Tailrace;
vector<Type> eEffectTurbineUnit4 = Tailrace;
vector<Type> eEffectTurbineUnit5 = Tailrace;
vector<Type> eEffectTurbineUnit6 = Tailrace;
vector<Type> eEffectSpill = Tailrace;
vector<Type> eEffectDischarge = Tailrace;
vector<Type> eReleaseEffect = Tailrace;
vector<Type> eOutput = Tailrace;
vector<Type> eRelease = Tailrace;
vector<Type> eTailrace = Tailrace;

Type nll = 0.0;

eTailrace(0) = Tailrace(0);

for(int i = 0; i < 6; i++){
  nll -= dnorm(bDepthUnit(i), Type(0), exp(sDepthUnit), true);

for(int i = 1; i < Tailrace.size(); i++){
  eTurbineUnit1(i) = Unit1(i) + 0.0001;
  eTurbineUnit2(i) = Unit2(i) + 0.0001;
  eTurbineUnit3(i) = Unit3(i) + 0.0001;
  eTurbineUnit4(i) = Unit4(i) + 0.0001;
  eTurbineUnit5(i) = Unit5(i) + 0.0001;
  eTurbineUnit6(i) = Unit6(i) + 0.0001;
  eDischarge(i) = eTurbineUnit1(i) + eTurbineUnit2(i) + eTurbineUnit3(i) + eTurbineUnit4(i) + eTurbineUnit5(i) + eTurbineUnit6(i) + Spill(i);
  eDepthTurbineUnit1(i) = bDepthUnit(0) + bDepthTurbineUnit(0) * eTurbineUnit1(i);
  eDepthTurbineUnit2(i) = bDepthUnit(1) + bDepthTurbineUnit(1) * eTurbineUnit2(i);
  eDepthTurbineUnit3(i) = bDepthUnit(2) + bDepthTurbineUnit(2) * eTurbineUnit3(i);
  eDepthTurbineUnit4(i) = bDepthUnit(3) + bDepthTurbineUnit(3) * eTurbineUnit4(i);
  eDepthTurbineUnit5(i) = bDepthUnit(4) + bDepthTurbineUnit(4) * eTurbineUnit5(i);
  eDepthTurbineUnit6(i) = bDepthUnit(5) + bDepthTurbineUnit(5) * eTurbineUnit6(i);
  eDepthSpill(i) = bDepthSpill0;
  eStratification(i) = (Forebay(i) - Intake(i)) / (707 - 690) ;
  eEffectTurbineUnit1(i) = eDepthTurbineUnit1(i) * eStratification(i);
  eEffectTurbineUnit2(i) = eDepthTurbineUnit2(i) * eStratification(i);
  eEffectTurbineUnit3(i) = eDepthTurbineUnit3(i) * eStratification(i);
  eEffectTurbineUnit4(i) = eDepthTurbineUnit4(i) * eStratification(i);
  eEffectTurbineUnit5(i) = eDepthTurbineUnit5(i) * eStratification(i);
  eEffectTurbineUnit6(i) = eDepthTurbineUnit6(i) * eStratification(i);
  eEffectSpill(i) = eDepthSpill(i) * eStratification(i);
  eEffectDischarge(i) = (eEffectTurbineUnit1(i) * eTurbineUnit1(i) + eEffectTurbineUnit2(i) * eTurbineUnit2(i) + eEffectTurbineUnit3(i) * eTurbineUnit3(i) + eEffectTurbineUnit4(i) * eTurbineUnit4(i) + eEffectTurbineUnit5(i) * eTurbineUnit5(i) + eEffectTurbineUnit6(i) * eTurbineUnit6(i) + eEffectSpill(i) * Spill(i)) / eDischarge(i);
  eReleaseEffect(i) = exp(bRelease + bReleaseAir * (Air(i) - Intake(i)));
  eOutput(i) = Intake(i) + eEffectDischarge(i);
  eRelease(i) = eOutput(i) + eReleaseEffect(i);
  eTailrace(i) = eRelease(i) + Consecutive(i) * invlogit(bRho) * (Tailrace(i-1) - eTailrace(i-1));
  nll -= dnorm(Tailrace(i), eTailrace(i), exp(sTailrace), true);
return nll;

Data Preparation

The fish and temperature data were provided by the Ktunaxa Nation. The discharge and elevation data were queried from the Columbia Basin Hydrological Database.

The data were cleaned and tidied using R version 4.0.1 (R Core Team 2015).

Length Cutoffs

Individuals were classified as fry (age-0), juvenile (age-1 and older subadults) or adult (sexually mature) based on the length cut-offs in Table 1.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on Bayesian estimation the reader is referred to McElreath (2016).

Unless indicated otherwise, the Bayesian analyses used normal and uniform prior distributions that were vague in the sense that they did not constrain the posteriors (Kery and Schaub 2011, 36). The posterior distributions were estimated from 1500 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of 3 chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that the potential scale reduction factor \(\hat{R} \leq 1.05\) (Kery and Schaub 2011, 40) and the effective sample size (Brooks et al. 2011) \(\textrm{ESS} \geq 150\) for each of the monitored parameters (Kery and Schaub 2011, 61).

The parameters are summarised in terms of the point estimate, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). The estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{mean}/\mathrm{sd}\) and the 95% CLs are the 2.5th and 97.5th percentiles. A p-value of 0.05 indicates that the lower or upper 95% CL is 0.

The results are displayed graphically by plotting the modeled relationships between particular variables and the response(s) with the remaining variables held constant. In general, continuous and discrete fixed variables are held constant at their mean and first level values, respectively, while random variables are held constant at their typical values (expected values of the underlying hyperdistributions) (Kery and Schaub 2011, 77–82). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% confidence/credible intervals (CIs, Bradford, Korman, and Higgins 2005).

The analyses were implemented using R version 4.0.1 (R Core Team 2018) and the mbr family of packages.

Model Descriptions

Body Condition

The annual variation in condition (body weight when accounting for body length) was estimated from the boat and backpack electrofishing captures using a mass-length model (He et al. 2008).

Key assumptions of the condition model include:

• Weight varies with body length as an allometric relationship, i.e., \(W = \alpha L^{\beta}\).
• \(\alpha\) varies with year.
• \(\beta\) varies with year.
• The residual variation in weight is log-normally distributed.

Preliminary analyses indicated that site and day of the year were not informative predictors of condition.

Relative Abundance

The annual variation in relative abundance was estimated from the boat count and catch data using an over-dispersed Poisson model (Kery and Schaub 2011). Lineal densities are by kilometre of river (as opposed to kilometre of bank).

Key assumptions of the relative abundance model include:

• Lineal density varies by period.
• Lineal density varies randomly with year.
• Lineal catch density is a fixed proportion of lineal count density.
• Expected counts (and catches) are the product of the count (catch) density and the length of river (half the length of bank) sampled.
• Observed counts (and catches) are described by a Poisson-gamma distribution.

Preliminary analyses indicated that site was not an informative predictor of lineal density.

The model does not distinguish between the abundance and observer efficiency, i.e., it estimates the count which is the product of the two. As such it is necessary to assume that changes in observer efficiency by year are negligible in order to interpret the estimates as relative abundance.

Water Temperature

Climatic variation can cause large differences in annual temperatures. Consequently, we explored the data for an effect of the additional turbines on the difference in the water temperature between the right versus left bank and when moving downstream. All apparent systematic differences were within the accuracy of the temperature loggers (\(\pm 0.2^{\circ}\text{C}\)).

Model Templates

.model{
  bWeightAlpha ~ dnorm(5, 2^-2)
  bWeightBeta ~ dnorm(3, 2^-2)
  bWeightAlphaYear[1] <- 0
  for(i in 2:nYear) {
    bWeightAlphaYear[i] ~ dnorm(0, 2^-2)
  }
  bWeightBetaYear[1] <- 0
  for(i in 2:nYear) {
    bWeightBetaYear[i] ~ dnorm(0, 2^-2)
  }
  sWeight ~ dnorm(0, 2^-2) T(0,)
  for (i in 1:length(Weight)) {
    eWeightAlpha[i] <- bWeightAlpha + bWeightAlphaYear[Year[i]]
    eWeightBeta[i] <- bWeightBeta + bWeightBetaYear[Year[i]]
    log(eWeight[i]) <- eWeightAlpha[i] + eWeightBeta[i] * Length[i]
    Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
  }

.model{
  bEfficiency <- 1
  bEfficiencyVisitType[1] <- 0
  for (i in 2:nVisitType) {
    bEfficiencyVisitType[i] ~ dunif(0, 1)
  }
  bDensity ~ dnorm(0, 5^-2)
  bDensityPeriod[1] <- 0
  for(i in 2:nPeriod) {
    bDensityPeriod[i] ~ dnorm(0, 2^-2)
  }
  sDensityYear ~ dnorm(0, 2^-2) T(0, )
  for(i in 1:nYear) {
    bDensityYear[i] ~ dnorm(0, sDensityYear^-2)
  }
  sDispersion ~ dnorm(0, 2^-2) T(0, )
  for (i in 1:length(Year)) {
    
    eEfficiency[i] <- bEfficiency - bEfficiencyVisitType[VisitType[i]]
    log(eDensity[i]) <- bDensity + bDensityPeriod[Period[i]] + bDensityYear[Year[i]]
    eAbundance[i] <- eDensity[i] * SiteLength[i] / 2
    eDispersion[i] ~ dgamma(1 / sDispersion^2, 1 / sDispersion^2)
    Count[i] ~ dpois(eAbundance[i] * eEfficiency[i] * eDispersion[i])
  }

Data Preparation

The 2017 survey data were provided by LGL Limited in the form of Excel tables. They were imported into the existing SQLite database. The data were clean and tidied (Wickham 2014) using R version 3.4.3 (R Core Team 2016). When there were repeated counts for a channel the higher of the two values was choosen. The mean daily discharge and water temperature values at Duncan River below Lardeau (DRL) were extracted from BC Hydro’s environmental database for the Kootenays, which is maintained by Poisson Consulting. The annual abundances of Kokanee spawning in Meadow Creek and the Lardeau River were provided by the Ministry of Forests, Lands and Natural Resource Operations.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The estimates were produced using JAGS (Plummer 2015). For additional information on Bayesian estimation the reader is referred to McElreath (2016).

Unless indicated otherwise, the Bayesian analyses used uninformative normal prior distributions (Kery and Schaub 2011, 36). The posterior distributions were estimated from 1500 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of 3 chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that \(\hat{R} \leq 1.1\) (Kery and Schaub 2011, 40) and \(\textrm{ESS} \geq 150\) for each of the monitored parameters (Kery and Schaub 2011, 61). Where \(\hat{R}\) is the potential scale reduction factor and \(\textrm{ESS}\) is the effective sample size.

The parameters are summarised in terms of the point estimate, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). The estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{mean}/\mathrm{sd}\) and the 95% CLs are the 2.5th and 97.5th percentiles. A p-value of 0.05 indicates that the lower or upper 95% CL is 0.

Where relevant, model adequacy was confirmed by examination of residual plots for the full model(s).

The results are displayed graphically by plotting the modeled relationships between particular variables and the response(s) with the remaining variables held constant. In general, continuous and discrete fixed variables are held constant at their mean and first level values, respectively, while random variables are held constant at their typical values (expected values of the underlying hyperdistributions) (Kery and Schaub 2011, 77–82). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% confidence/credible intervals (CIs, Bradford, Korman, and Higgins 2005).

The analyses were implemented using R version 3.4.3 (R Core Team 2015) and the mbr family of packages.

Model Descriptions

Model Templates

Observer Efficiency

model{

  bEfficiency ~ dunif(0, 3)
  sDispersion ~ dunif(0, 5)

  for (i in 1:nObs){
    eAerial[i] <- Ground[i] * bEfficiency
    eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
    Aerial[i] ~ dpois(eAerial[i] * eDispersion[i])
  }
..

Template 1. Model definition.

Area-Under-The-Curve

model{
  bAbundance ~ dnorm(10, 5^-2)
  bPeakTiming ~ dunif(-14, 14)
  bDuration ~ dnorm(2, 1^-2)

  bEfficiency ~ dnorm(1.40, 0.38^-2) T(0.89, 2.40)
  bResidenceTime ~ dunif(7, 14)

  sAbundanceYear ~ dunif(0, 5)
  sPeakTiming ~ dunif(0, 14)
  for(i in 1:nYear){
    bAbundanceYear[i] ~ dnorm(0, sAbundanceYear^-2)
    bPeakTimingYear[i] ~ dnorm(0, sPeakTiming^-2)
  }

  bSD ~ dnorm(0, 10^-2)
  bSDAbundance ~ dnorm(0, 5^-2)
  for(i in 1:length(Count)){

    log(eAbundance[i]) <- bAbundance + bAbundanceYear[Year[i]]
    ePeakTiming[i] <- bPeakTiming + bPeakTimingYear[Year[i]]
    log(eDuration[i]) <- bDuration

    eArrived[i] <- pnorm(Dayte[i], (ePeakTiming[i] - bResidenceTime/2), eDuration[i]^-2)
    eDeparted[i] <- pnorm(Dayte[i], (ePeakTiming[i] + bResidenceTime/2), eDuration[i]^-2)
    eSpawners[i] <- (eArrived[i] - eDeparted[i]) * eAbundance[i]

    eEfficiency[i] <- bEfficiency

    log(eSD[i]) <- bSD + bSDAbundance * eAbundance[i]

    Count[i] ~ dnorm(eSpawners[i] * eEfficiency[i], eSD[i]^-2)
  }
..

Template 2. Model definition.