Growth

Lower Columbia River Fish Population Indexing 2023

The suggested citation for this analytic appendix is:

Thorley, J.L. (2024) Lower Columbia River Fish Population Indexing 2023. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/25988200.

Background

In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.

The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:

What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.

Lower Columbia River Fish Population Indexing 2022

The suggested citation for this analytic appendix is:

Thorley, J.L., Dalgarno, S., Lyons, S. (2023) Lower Columbia River Fish Population Indexing 2022. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1392245070.

Background

In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.

The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:

What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.

Lower Columbia River Fish Population Indexing 2021

The suggested citation for this analytic report is:

Thorley, J.L. & Hussein, N. (2022) Lower Columbia River Fish Population Indexing 2021. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1314086833.

Background

In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.

The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:

What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River? What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River? The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.

Horse Lake Exploitation Analysis 2021

The suggested citation for this analytic appendix is:

Thorley, J.L. & Pearson, A.D. (2022) Horse Lake Exploitation Analysis 2021. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/556371688.

Background

Horse Lake supports a recreational fishery for Lake Trout (Salvelinus namaycush). To estimate their natural and fishing mortality, individuals were caught by angling and tagged with acoustic transmitters and/or external reward tags. Acoustically tagged fish were detected by a series of seven receivers distributed through the lake as well as annual or biannual mobile detection surveys. The mobile detection survey was conducted by deploying receivers for approximately 10 minutes at 90 locations arranged on a 400 m x 400 m grid. The external reward tags included a phone number for anglers to report their recaptures. In 2009 and 2021, a summer profundal index netting (SPIN) survey was done on the lake.

All the data were provided by the Ministry of Forests, Lands and Natural Resource Operations (MFLNRO) Cariboo Region.

Data Preparation

The data were prepared for analysis using R version 4.2.0 (R Core Team 2021). Detections were aggregated daily, where for each transmitter the receiver with the most detections was chosen.

The Seasons were Winter (December - February), Spring (March - May), Summer (June - August) and Autumn (September - November).

For each mobile survey each detected fish was assigned a location corresponding to the centroid of its detection locations weighted by the detection rates. The mobile receiver deployment times were corrected by increasing or decreasing the deployment times for each receiver on each day by up to three minutes to maximize the number of detections within deployments. Mobile surveys conducted within 3 months of each other were combined into a single survey.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The estimates were produced using JAGS (Plummer 2003) and STAN (Carpenter et al. 2017). 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 (Chow and Greenland 2019) is the Shannon transform (-log to base 2) of the corresponding p-value (Kery and Schaub 2011; Greenland and Poole 2013). A surprisal value of 4.3 bits, which is equivalent to a p-value of 0.05 indicates that the surprise would be equivalent to throwing 4.3 heads in a row. The condition that non-essential explanatory variables have s-values \(\geq\) 4.3 bits provides a useful model selection heuristic (Kery and Schaub 2011).

Model adequacy was assessed via posterior predictive checks (Kery and Schaub 2011). More specifically, the first four central moments (mean, variance, skewness and kurtosis) for the deviance residuals were compared to the expected values by simulating new residuals. In this context the s-value indicates how surprising each metric is given the estimated posterior probability distribution for the residual variation.

Where computationally practical, the sensitivity of the parameters to the choice of prior distributions was evaluated by increasing the standard deviations of all normal, half-normal and log-normal priors by an order of magnitude and then using \(\hat{R}\) to test whether the samples where drawn from the same posterior distribution (Thorley and Andrusak 2017).

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.2.0 (R Core Team 2021) and the mbr family of packages.

Model Descriptions

Condition

The expected weight, \(W\), of a fish of a given length, \(L\), was estimated from the data using an allometric mass-length model (He et al. 2008)

\[ W = \alpha L^\beta\]

Key assumptions of the condition model include:

  • The prior distribution on the power term (\(\beta\)) is a normal distribution with a mean of 3.18 and a standard deviation of 0.19 (McDermid, Shuter, and Lester 2010).
  • The weight varies by day of the year.
  • The residual variation in weight is log-normally distributed.

Preliminary analysis indicated year as a random effect and capture method (angling vs SPIN) were not important predictors of condition.

Growth and Spawning

The growth curve and length at spawning for each ecotype were estimated from the length, scale age, and spawning data using an integrated model. The model consists of a biphasic Von Bertalanffy (von Bertalanffy 1938) growth curve and length at spawning model.

The Von Bertalanffy growth curve (von Bertalanffy 1938) for the small ecotype is

\[ L = L_{\infty_{1}} \cdot (1 - \exp(-k_1 \cdot A ))\]

where \(A\) is the scale age in years.

The large ecotype transitions to the second growth phase at the transitions age (\(A_T\)), where \(L_T\) is the transition length (mm)

\[ L_T = L_{\infty_{1}} \cdot (1 - \exp(-k_1 \cdot A_T )) \]

The Von Bertalanffy growth curve (von Bertalanffy 1938) in the second growth phase the large ecotype

\[ L = L_T + (L_{\infty_{2}} - L_T) \cdot (1 - \exp(-k_2 \cdot (A - A_T)))\]

The probability of spawning at length \(L\) is determined by

\[S = \frac{L^{S_p}}{L_s^{S_p} + L^{S_p}} \cdot es\]

Key assumptions of the model include:

  • \(Ls\) varies by ecotype
  • \(Es\) varies by ecotype
  • Transition age, \(A_T\), varies randomly by fish
  • The residual variation in length is log normally distributed.

The sensitivity analysis indicated that the estimate of bBigEs is sensitive to the choice of prior. Increasing the standard deviation of the prior by an order of magnitude increased the estimate from 1.71 (95% CI 0.06 - 4.15) to 8.97 (95% CI 0.153 - 42.0).

Mobile

A fish was classified as being in the lake at the time of the survey if it was detected 1 or more times. A fish was classified as having moved if the weighted centriod had shifted more then 800 m since the previous survey.

Survival

The natural mortality and capture probability for individuals \(\geq\) 300 mm was estimated using a Bayesian individual state-space survival model (Thorley and Andrusak 2017) with three month intervals. The survival model incorporated natural mortality, fixed and mobile acoustic detection and movement, T-bar tag loss, recapture and reporting. Key assumptions of the survival model include:

  • The natural mortality varies by season.
  • The fishing mortality varies by year.

The survival model treats all recaptured fish as if they had been retained. Only individuals with a fork length (FL) \(\geq\) 300 mm that were tagged with an acoustic tag and/or $100 and $10 reward tags were included in the survival analysis.

Preliminary analysis indicated that there was not a clear difference between the natural mortality and fishing mortality rates for individuals \(\geq\) 500 mm.

Yield-Per-Recruit

The optimal capture rate (to maximize number of individuals harvested assuming 100% retention) was calculated using yield-per-recruit analysis (Walters and Martell 2004). A yield-per-recruit analysis was performed for the two different Lake Trout ecotypes separately and in combination. The \(Rk\) of 25 from Myers et al. (Myers, Bowen, and Barrowman 1999) was converted into the equivalent egg to age-1 survival of 0.025 based on the large ecotype life history

Model Templates

Condition

.model {
  bWeight ~ dnorm(-12, 2^-2) 
  bWeightLength ~ dnorm(3.18, 0.19^-2)
  bDayte ~ dnorm(0, 2^-2)
  bDayte2 ~ dnorm(0, 2^-2)
  sWeight ~ dnorm(0, 1^-2) T(0,)
  for(i in 1:length(LogLength)) {
    eWeightLength[i] <- bWeightLength + bDayte * Dayte[i] + bDayte2 * Dayte[i]^2 
    eWeight[i] <- bWeight + eWeightLength[i] * LogLength[i] 
    Weight[i] ~ dlnorm(eWeight[i], sWeight^-2)
  }

Block 1. Condition model description.

Growth and Spawning

.model {
  bT0 <- 0
               
  bLInf1 ~ dnorm(6, 2^-2)
  bLInf2  ~ dnorm(7, 2^-2)
  bK1 ~ dnorm(0, 2^-2)
  bK2 ~ dnorm(0, 2^-2)             
  bTransitionAge ~ dnorm(2, 1^-2)
  sTransitionAgeFish ~ dnorm(0, 2^-2) T(0,)
  bSmall ~ dnorm(0, 2^-2) 
  sLength ~ dnorm(0, 2^-2) T(0, )
   
  Ls ~ dnorm(6, 2^-2)
  Sp ~ dnorm(20, 20^-2) T(0,)
  Es ~ dnorm(0, 2^-2)
  bBigLs ~ dnorm(0, 2^-2)
  bBigEs ~ dnorm(0, 2^-2)
  for (i in 1:nFishID) {
    bTransitionAgeFish[i] ~ dnorm(-sTransitionAgeFish^2/2, sTransitionAgeFish^-2)
  }
  for (i in 1:nObs) {
    logit(eSmall[i]) <- bSmall
    log(eLInf1[i]) <- bLInf1 
    log(eLInf2[i]) <- bLInf2
    log(eK1[i]) <- bK1
    log(eK2[i]) <- bK2  
    log(eLs[i]) <- Ls  + bBigLs * (1 - SmallEcotype[i])
    logit(eEs[i]) <- Es + bBigEs * (1 - SmallEcotype[i])
    eSpawning[i] <- ( (Length[i] ^ Sp) / (eLs[i]^Sp + Length[i]^Sp) ) * eEs[i]
    eSmallLength[i] <- eLInf1[i] * (1 - exp(-eK1[i] * (Age[i] - bT0))) 
    
    log(eTransitionAge[i]) <- bTransitionAge + bTransitionAgeFish[FishID[i]]
    eTransitionLength[i] <- eLInf1[i] * (1 - exp(-eK1[i] * (eTransitionAge[i] - bT0))) 
    eBigLength[i] <- eSmallLength[i] * step(eTransitionAge[i] - Age[i]) + step(Age[i] - eTransitionAge[i] - 0.01) * (eTransitionLength[i] + (eLInf2[i] - eTransitionLength[i]) * (1 - exp(-eK2[i] * (Age[i] - eTransitionAge[i]))))
    eLength[i] <-  eSmallLength[i] * SmallEcotype[i] + eBigLength[i] * (1 - SmallEcotype[i])
    Spawning[i] ~ dbin(eSpawning[i], 1)
    SmallEcotype[i] ~ dbin(eSmall[i], 1)
    Length[i] ~ dlnorm(log(eLength[i]), sLength^-2)
  }

Block 2. Growth model description.

Survival

.model{
  bMortality ~ dnorm(-3, 3^-2)
  bTagLoss ~ dbeta(1, 1)
  bMoved ~ dbeta(1, 1)
  bMobileDetected ~ dbeta(1, 1)
  bMobileMoved ~ dbeta(1, 1)
  bRecaptured ~ dnorm(0, 2^-2)
  bReported ~ dbeta(1, 1)
  bDetectedAlive ~ dbeta(1, 1)
  bDieNear ~ dbeta(1, 1)
  bDetectedDeadNear ~ dbeta(10, 1)
  bDetectedDeadFar ~ dbeta(1, 10)
  bMortalitySeason[1] <- 0
  for (i in 2:nSeason) {
    bMortalitySeason[i] ~ dnorm(0, 2^-2)
  }
  sRecapturedAnnual ~ dnorm(0, 2^-2) T(0,)
  for(i in 1:nAnnual) {
    bRecapturedAnnual[i] ~ dnorm(0, sRecapturedAnnual^-2)
  }

  for (i in 1:nCapture){
    eDetectedAlive[i] <- bDetectedAlive
    eDieNear[i] ~ dbern(bDieNear)
    eDetectedDeadNear[i] <- bDetectedDeadNear
    eDetectedDeadFar[i] <- bDetectedDeadFar
    eDetectedDead[i] <- eDieNear[i] * eDetectedDeadNear[i] + (1 - eDieNear[i]) * eDetectedDeadFar[i]
    logit(eMortality[i,PeriodCapture[i]]) <- bMortality + bMortalitySeason[Season[PeriodCapture[i]]]
    eTagLoss[i,PeriodCapture[i]] <- bTagLoss
    eDetected[i,PeriodCapture[i]] <- Alive[i,PeriodCapture[i]] * eDetectedAlive[i] + (1-Alive[i,PeriodCapture[i]]) * eDetectedDead[i]
    eMoved[i,PeriodCapture[i]] <- bMoved
    eMobileDetected[i,PeriodCapture[i]] <- bMobileDetected
    eMobileMoved[i,PeriodCapture[i]] <- bMobileMoved
    
    logit(eRecaptured[i,PeriodCapture[i]]) <- bRecaptured + bRecapturedAnnual[Annual[PeriodCapture[i]]]
    eReported[i,PeriodCapture[i]] <- bReported
    InLake[i,PeriodCapture[i]] <- 1
    Alive[i,PeriodCapture[i]] ~ dbern(1-eMortality[i,PeriodCapture[i]])
    TBarTag100[i,PeriodCapture[i]] ~ dbern(1-eTagLoss[i,PeriodCapture[i]])
    TBarTag10[i,PeriodCapture[i]] ~ dbern(1-eTagLoss[i,PeriodCapture[i]])
    Detected[i,PeriodCapture[i]] ~ dbern(Monitored[i,PeriodCapture[i]] * eDetected[i,PeriodCapture[i]])
    Moved[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Monitored[i,PeriodCapture[i]] * eMoved[i,PeriodCapture[i]])
    MobileDetected[i,PeriodCapture[i]] ~ dbern(Monitored[i,PeriodCapture[i]] * eMobileDetected[i,PeriodCapture[i]] * MobileSurvey[i,PeriodCapture[i]])
    MobileMoved[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * MobileDetected[i,PeriodCapture[i]] * eMobileMoved[i,PeriodCapture[i]])  
    Recaptured[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * eRecaptured[i,PeriodCapture[i]])
    Reported[i,PeriodCapture[i]] ~ dbern(Recaptured[i,PeriodCapture[i]] * eReported[i,PeriodCapture[i]] * step(TBarTag100[i,PeriodCapture[i]] + TBarTag10[i,PeriodCapture[i]] - 1))
  for(j in (PeriodCapture[i]+1):nPeriod) {
    logit(eMortality[i,j]) <- bMortality + bMortalitySeason[Season[j]]
    eTagLoss[i,j] <- bTagLoss

    eDetected[i,j] <- Alive[i,j] * eDetectedAlive[i] + (1-Alive[i,j]) * eDetectedDead[i]
    eMoved[i,j] <- bMoved
    eMobileDetected[i,j] <- bMobileDetected
    eMobileMoved[i,j] <- bMobileMoved
    logit(eRecaptured[i,j]) <- bRecaptured + bRecapturedAnnual[Annual[j]]
    eReported[i,j] <- bReported
    
    InLake[i,j] ~ dbern(InLake[i,j-1] * (1-Recaptured[i,j-1]))
    Alive[i,j] ~ dbern(Alive[i,j-1] * (1-Recaptured[i,j-1]) * (1-eMortality[i,j-1]))
    TBarTag100[i,j] ~ dbern(TBarTag100[i,j-1] * (1-Recaptured[i,j-1]) * (1-eTagLoss[i,j]))
    TBarTag10[i,j] ~ dbern(TBarTag10[i,j-1] * (1-Recaptured[i,j-1]) * (1-eTagLoss[i,j]))
    Detected[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eDetected[i,j])
    Moved[i,j] ~ dbern(Alive[i,j] * Monitored[i,j] * eMoved[i,j])
    MobileDetected[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eMobileDetected[i,j] * MobileSurvey[i,j])
    MobileMoved[i,j] ~ dbern(Alive[i,j] * MobileDetected[i,j] * eMobileMoved[i,j])
    Recaptured[i,j] ~ dbern(Alive[i,j] * eRecaptured[i,j])
    Reported[i,j] ~ dbern(Recaptured[i,j] * eReported[i,j] * step(TBarTag100[i,j] + TBarTag10[i,j] - 1))}}
  }

Block 3. Survival model description.

Lower Columbia River Fish Population Indexing 2020

The suggested citation for this analytic report is:

Thorley, J.L. & Hussein, N. (2021) Lower Columbia River Fish Population Indexing 2020. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/392669554.

Background

In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.

The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:

What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.

Mica Dam Expansion Water Temperature and Fish Indexing Study 2020

The suggested citation for this analytic appendix is:

Thorley, J.L. & Amies-Galonski, E. (2021) Mica Dam Expansion Water Temperature and Fish Indexing Study 2020. A Poisson Consulting Analysis Appendix. URL: https://www.poissonconsulting.ca/f/955630057.

Background

The Mica Tailrace Fish Indexing Study is a multi-year program to estimate the effects of the addition of two new turbines (Mica 5 and 6) on the ichyofauna and thermal regime in the 2.5 km of the Columbia River downstream of Mica Dam. A single year of fish indexing data (2008) was also available from a previous program. As per the Terms of Reference (TOR) the relative abundance, condition and spatial distribution of the fish populations was assessed.

Mica 5 became operational on January 28th 2015 and Mica 6 became operational on December 22nd 2015.

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

Tailrace

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}\)).

Forebay and Tailrace

In 2020 the BC Hydro thermistor data for the Mica Dam forebay were analysed using Maximum Likelihood (Millar 2011) and TMB (Kristensen et al. 2016) to determine the extent to which discharge from the turbines and the air temperature influences the water temperature in the tailrace.

Key assumptions of the hourly temperature model include:

  • The effective depth of each unit depends on the unit (as a random effect) and its discharge (as a fixed effect).
  • The water temperature effect of each unit depends on its effective depth and the stratification in the forebay.
  • The water temperature is affected by release from the dam and the difference in air temperature and the water temperature.
  • The temporal autocorrelation is described a first order moving average process.
  • The residual water temperature at the tailrace temperature logger is normally distributed.

Despite the incorporation of a first order moving average process the residual variation in the hourly tailrace water temperature was strongly autocorrelated. Consequently CIs are not provided for the coefficients or estimates.

Model Templates

Condition

.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)
  }

Block 1.

Relative Abundance

.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])
  }

Block 2.

Forebay and Tailrace

.#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;

Block 3. Model description.

Horse Lake Exploitation Analysis 2020

The suggested citation for this analytic appendix is:

Thorley, J.L., Dalgarno, S. & Pearson, A. (2021) Horse Lake Exploitation Analysis 2020. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1065913173.

Background

Horse Lake supports a recreational fishery for Lake Trout (Salvelinus namaycush). To estimate their natural and fishing mortality, individuals were caught by angling and tagged with acoustic transmitters and/or external reward tags. Acoustically tagged fish were detected by a series of seven receivers distributed through the lake as well as annual mobile receiver surveys. The external reward tags included a phone number for anglers to report their recaptures. In 2019, eight Rainbow Trout with acoustic transmitters and external reward tags were transplanted from Greenlee Lake.

All the data were provided by the Ministry of Forests, Lands and Natural Resource Operations (MFLNRO) Cariboo Region.

Lower Columbia River Fish Population Indexing 2019

The suggested citation for this analytic report is:

Thorley, J.L. (2020) Lower Columbia River Fish Population Indexing 2019. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1566579922.

Background

In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.

The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:

What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?

The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.

Mica Dam Expansion Water Temperature and Fish Indexing Study 2019

The suggested citation for this analytic appendix is:

Thorley, J.L. & Amies-Galonski, E. (2020) Mica Dam Expansion Water Temperature and Fish Indexing Study 2019. A Poisson Consulting Analysis Appendix. URL: http://www.poissonconsulting.ca/f/1134568871.

Background

The Mica Tailrace Fish Indexing Study is a multi-year program to estimate the effects of the addition of two new turbines (Mica 5 and 6) on the ichyofauna and thermal regime in the 2.5 km of the Columbia River downstream of Mica Dam. A single year of fish indexing data (2008) was also available from a previous program. As per the Terms of Reference (TOR) the relative abundance, condition and spatial distribution of the fish populations was assessed. In addition, changes in the species evenness were also estimated.

Mica 5 became operational on January 28th 2015 and Mica 6 became operational on December 22nd 2015.

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

Condition

.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)
  }

Block 1.

Relative Abundance

.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])
  }

Block 2.

Middle Columbia River Fish Indexing Analysis 2019

The suggested citation for this analytic report is:

Thorley, J.L. and Amies-Galonski E. (2020) Middle Columbia River Fish Indexing Analysis 2019. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1050384286.

Background

The key management questions to be addressed by the analyses are:

  1. Is there a change in abundance of adult life stages of fish using the Middle Columbia River (MCR) that corresponds with the implementation of a year-round minimum flow?
  2. Is there a change in growth rate of adult life stages of the most common fish species using the MCR that corresponds with the implementation of a year-round minimum flow?
  3. Is there a change in body condition (measured as a function of relative weight to length) of adult life stages of fish using the MCR that corresponds with the implementation of a year-round minimum flow?
  4. Is there a change in spatial distribution of adult life stages of fish using the MCR that corresponds with the implementation of a year-round minimum flow?

Other objectives include the estimation of species richness and diversity. The year-round minimum flow was implemented in the winter of 2010 at the same time that a fifth turbine was added.