Lower Columbia River Rainbow Trout Spawning Analysis 2014

The main intepretive report by Mountain Waters Research and Poisson Consulting, which was prepared for BC Hydro, is available from BC Hydro.

The suggested citation for this online appendix is:

Thorley, J.L. (2015) Lower Columbia River Rainbow Trout Spawning Analysis 2014. URL: http://www.poissonconsulting.ca/f/1157179155.

Background

Each spring in the Lower Columbia River (LCR) below Hugh L. Keenleyside Dam (HLK) and in the Lower Kootenay River (LKR) below Brilliant Dam, thousands of Rainbow Trout spawn. Since 1992, BC Hydro has stabilized the spring discharge releases from HLK to protect Rainbow Trout redds from dewatering.

The primary goal of the current analysis to answer the following management question:

Does the implementation of Rainbow Trout Spawning Protection Flows over the course of the monitoring period lead to an increase in the relative abundance of Rainbow Trout spawning in the LCR downstream of HLK dam?

Methods

Data Preparation

The Rainbow Trout fish and redd aerial count data for the LCR and LKR were collected by Mountain Water Research and databased by Gary Pavan. The age-1 Rainbow Trout abundance data were provided by the Fish Population Indexing Program which is led by Golder Associates.

The study area was divided into three sections: the LCR above the LKR, the LKR and the LCR below the LKR. Redd and spawner counts upstream of Norns Creek Fan and downstream of Genelle were excluded from the section totals because they constitute less than 0.1% of the total count and were not surveyed in all years. The redd and spawner counts for the Right Upstream Bank above Robson Bridge were also excluded as they appear to be primarily driven by viewing conditions (and consitute less than 2.5% of the total). Simultaneous declines in both spawners and redds for a particular section were inferred to be caused by poor viewing conditions (turbidity) and the affected spawner and redd section counts were excluded from any subsequent analyses.

Statistical Analysis

Hierarchical Bayesian models were fitted to the count data using R version 3.1.1 (Team, 2013) and JAGS 3.4.0 (Plummer, 2012) which interfaced with each other via jaggernaut 1.8.2 . For additional information on hierarchical Bayesian modelling in the BUGS language, of which JAGS uses a dialect, the reader is referred to Kery and Schaub (2011) pages 41-44.

Unless specified, the models assumed vague (low information) prior distributions . The posterior distributions were estimated from a minimum of 1,000 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of three chains . Model convergence was confirmed by ensuring that Rhat was less than 1.1 for each of the parameters in the model . Model adequacy was confirmed by examination of residual plots.

The posterior distributions of the fixed parameters are summarised in terms of a point estimate (mean), lower and upper 95% credible limits (2.5th and 97.5th percentiles), the standard deviation (SD), percent relative error (half the 95% credible interval as a percent of the point estimate) and significance .

Model selection was achieved by dropping insignificant explanatory parameters and uninformative hyperparameters. For the purposes of the current analysis an explanatory parameter was considered to be insignificant if its significance was \(\geq\) 0.05 and a hyperparameter was considered to be uninformative if its percent relative error was \(\geq\) 80%.

The results are displayed graphically by plotting the modeled relationships between particular variables and the response with 95% credible intervals (CRIs) 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) . Where informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% CRIs (Bradford, Korman, and Higgins, 2005). Plots were produced using the ggplot2 R package .

Area-Under-the-Curve

The spawner abundance and spawn timing was estimated from the aerial fish and redd counts for the three sections using an Area-Under-the-Curve (AUC) model.

Key assumptions of the AUC model include:

The expected emergence timing was calculated from the estimated spawn timing and the surface water temperature under the assumption that Rainbow Trout embryos emerge after 480 accumulated thermal units (ATUs).

The proportions of spawners at each site were used to calculate the Shannon Index, an information-theoretic measure of the diversity in the abundance distribution of a resource (Krebs, 1999). In the current context, the Shannon Index takes into account both the number of spawning sites and how the spawning activity is distributed among these, with a higher index indicating a greater spatial distribution of spawning.

The Shannon Index (\(H\)) is given by \[ H = - \sum{p_i log(p_i)}\] where \(p_i\) is the proportion of the spawning activity at the ith location.

Stock-Recruitment

The relationship between the number of spawners and the number of age-1 recruits in the fall of the following year was estimated using a Beverton-Holt stock-recruitment model.

Key assumptions of the stock-recruitment model include:

The prior probability distribution mean of 90 for R0 was based on an average of 2,900 eggs per female spawner, a 50:50 sex ratio, 50% egg survival, 50% post-emergence fall survival, 50% overwintering survival and 50% summer survival.

Model Code

The JAGS model code, which uses a series of naming conventions, is presented below.

Area-Under-The-Curve

Variable/Parameter Description
bPeakSpawnerArrivalTiming Timing of peak of spawner arrival
bReddObserverEfficiency Redd observer efficiency
bReddPerSpawner Number of redds per spawner
bReddResidenceTime Redd residence time
bSpawnerAbundanceSite[i] Intercept for spawner abundance at ith site
bSpawnerObserverEfficiency Spawner observer efficiency
bSpawnerResidenceTime Spawner residence time
Dayte[i] Day of the year of the ith count
Fish[i] Number of fish observed on ith count
Redds[i] Number of redds observed on ith count
sFishDispersion SD of overdispersion for Fish
Site[i] Site of the ith count
sPeakSpawnerArrivalTimingYear SD of effect of Year on bPeakSpawnerArrivalTiming
sReddsDispersion SD of overdispersion for Redds
sSpawnDuration SD of the duration of spawning
sSpawnerAbundanceSiteYear SD of effect of Site within Year on bSpawnerAbundanceSite
sSpawnerAbundanceYear SD of effect of Year on bSpawnerAbundanceSite
Year[i] Year of the ith count
Area-Under-The-Curve - Model1
model {

  sSpawnDuration ~ dunif(14, 42)

  bPeakSpawnerArrivalTiming ~ dnorm(130, 14^-2)

  bSpawnerResidenceTime ~ dunif(14, 21)
  bSpawnerObserverEfficiency ~ dunif(0.9, 1.1)

  bReddObserverEfficiency ~ dunif(0.9, 1.1)
  bReddResidenceTime <- 35

  bReddPerSpawner ~ dunif(0, 4)

  sPeakSpawnerArrivalTimingYear ~ dunif(0, 28)
  sSpawnerAbundanceYear ~ dunif(0, 2)
  for (i in 1:nYear) {
    bPeakSpawnerArrivalTimingYear[i] ~ dnorm (0, sPeakSpawnerArrivalTimingYear^-2)
    bSpawnerAbundanceYear[i] ~ dnorm (0, sSpawnerAbundanceYear^-2)
  }

  sSpawnerAbundanceSiteYear ~ dunif (0, 2)
  for (i in 1:nSite) {
    bSpawnerAbundanceSite[i] ~ dnorm (5, 5^-2)
    for (j in 1:nYear) {
      bSpawnerAbundanceSiteYear[i, j] ~ dnorm (0, sSpawnerAbundanceSiteYear^-2)
    }
  }

  sFishDispersion ~ dunif(0, 2)
  sReddDispersion ~ dunif(0, 2)
  for (i in 1:length(Fish)) {
    ePeakSpawnerArrivalTiming[i] <- bPeakSpawnerArrivalTiming + bPeakSpawnerArrivalTimingYear[Year[i]]

    log(eSpawnerAbundance[i]) <- bSpawnerAbundanceSite[Site[i]] + bSpawnerAbundanceYear[Year[i]] + bSpawnerAbundanceSiteYear[Site[i], Year[i]]

    eFish[i] <- (pnorm(Dayte[i], ePeakSpawnerArrivalTiming[i], sSpawnDuration^-2) - pnorm(Dayte[i], ePeakSpawnerArrivalTiming[i] + bSpawnerResidenceTime, sSpawnDuration^-2)) * eSpawnerAbundance[i] * bSpawnerObserverEfficiency

    eRedds[i] <- (pnorm(Dayte[i], ePeakSpawnerArrivalTiming[i], sSpawnDuration^-2) - pnorm(Dayte[i], ePeakSpawnerArrivalTiming[i] + bReddResidenceTime, sSpawnDuration^-2)) * eSpawnerAbundance[i] * bReddPerSpawner * bReddObserverEfficiency

    eFishDispersion[i] ~ dgamma(1 / sFishDispersion^2, 1 / sFishDispersion^2)
    eReddDispersion[i] ~ dgamma(1 / sReddDispersion^2, 1 / sReddDispersion^2)

    Fish[i] ~ dpois(eFish[i] *  eFishDispersion[i])
    Redds[i] ~ dpois(eRedds[i] *  eReddDispersion[i])
  }
} 

Stock-Recruitment

Variable/Parameter Description
eRecruits Expected number of Recruits
k Maximum Recruits
R0 Maximum Recruits per Stock
Recruits Number of age-1 individuals the following year
sRecruits SD of log-normally distributed residual variation in Recruits
Stock Number of spawners
Stock-Recruitment - Model1
model {

  R0 ~ dnorm(2900 * 0.5^5, 50^-2) T(0, )

  k ~ dlnorm(10, 5^-2)

  sRecruits ~ dunif(0, 5)  
  for (i in 1:length(Stock)) {
    eRecruits[i] <- R0 * Stock[i] / (1 + Stock[i] * (R0 - 1) / k)
    Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
  }
}

Results

Model Parameters

The posterior distributions for the fixed (Kery and Schaub 2011 p. 75) parameters in each model are summarised below.

Area-Under-The-Curve

Parameter Estimate Lower Upper SD Error Significance
bPeakSpawnerArrivalTiming 126.1914 120.9502 131.8695 2.7488 4 0
bReddObserverEfficiency 0.9955 0.9046 1.0934 0.0561 9 0
bReddPerSpawner 0.7661 0.6157 0.9294 0.0851 20 0
bReddResidenceTime 35.0000 35.0000 35.0000 0.0000 0 0
bSpawnerAbundanceSite[1] 7.3473 6.9451 7.7525 0.2086 5 0
bSpawnerAbundanceSite[2] 6.6450 6.2269 7.0652 0.2122 6 0
bSpawnerAbundanceSite[3] 7.8584 7.4406 8.2831 0.2117 5 0
bSpawnerObserverEfficiency 1.0061 0.9060 1.0948 0.0564 9 0
bSpawnerResidenceTime 17.0225 14.1652 20.5962 1.8557 19 0
sFishDispersion 0.7832 0.7241 0.8411 0.0303 7 0
sPeakSpawnerArrivalTimingYear 8.2512 5.3229 13.1978 2.0025 48 0
sReddDispersion 0.2851 0.2589 0.3135 0.0144 10 0
sSpawnDuration 27.8910 26.6355 29.2164 0.6708 5 0
sSpawnerAbundanceSiteYear 0.2206 0.1575 0.3013 0.0373 33 0
sSpawnerAbundanceYear 0.6642 0.4296 1.0110 0.1593 44 0
Rhat Iterations
1.02 8e+05

Stock-Recruitment

Parameter Estimate Lower Upper SD Error Significance
k 2.350e+04 1.971e+04 2.750e+04 2009.0000 17 0
R0 1.084e+02 3.689e+01 2.012e+02 43.1080 76 0
sRecruits 3.208e-01 2.157e-01 4.815e-01 0.0706 41 0
Rhat Iterations
1.09 10000

Figures

Discharge

figures/discharge/pre-post.png
Figure 1. Mean, minimum and maximum hourly discharge by date, location and pre (1999-2006) versus post (2007-2014) period.

Temperature

figures/temperature/temperature.png
Figure 2. Mean, minimum and maximum daily water temperature by date, location and pre (1999-2006) versus post (2007-2014) period. The time series were incomplete.

Acoustic Detections

figures/detections/fan-rainbow.png
Figure 3. Rainbow Trout detections at Norns Creek fan by date, year and diel period. The grey area indicates the period of receiver deployment.
figures/detections/creek-rainbow.png
Figure 4. Rainbow trout detections in Norns Creek by date and location.
figures/detections/residence.png
Figure 5. Rainbow Trout residence times at Norns Creek fan by year where residence is considered any detection within a 24 hr period.

Area-Under-The-Curve

figures/auc/upstream.png
Figure 6. Aerial fish and redds counts and predicted values in the LCR above the LKR by date and year.
figures/auc/kootenay.png
Figure 7. Aerial fish and redds counts and predicted values in the LKR by date and year.
figures/auc/downstream.png
Figure 8. Aerial fish and redds counts and predicted values in the LCR below the LKR by date and year.
figures/auc/year.png
Figure 9. Estimated total spawner abundance by year (with 95% CRIs).
figures/auc/spawn_timing.png
Figure 10. Estimated start, peak and end spawn timing by year (with 95% CRIs).
figures/auc/emergence_timing.png
Figure 11. Estimated start, peak and end emergence timing by year (with 95% CRIs).
figures/auc/peak-percent.png
Figure 12. Percent of peak spawner count by river kilometre and year.
figures/auc/shannon.png
Figure 13. Shannon index for the spatial distribution of spawners by year.
figures/auc/year_redds.png
Figure 14. Estimated percent of redds dewatered by year (with 95% CRIs).

Stock-Recruitment

figures/sr/subadult.png
Figure 15. Estimated age-1 recruits by year (with 95% CRIs) from LCR Fish Population Indexing.
figures/sr/sr.png
Figure 16. Estimated stock-recruitment curve (with 95% CRIs).

Norns Creek

figures/nornscreek/nornsabundance.png
Figure 17. Estimated spawner abundance in Norns Creek by year.
figures/nornscreek/nornsdistribution.png
Figure 18. Distribution of spawning fish in Norns Creek by year.

References

[1] M. J. Bradford, J. Korman and P. S. Higgins. “Using confidence intervals to estimate the response of salmon populations (Oncorhynchus spp.) to experimental habitat alterations”. In: Canadian Journal of Fisheries and Aquatic Sciences 62.12 (2005), pp. 2716-2726. ISSN: 0706-652X, 1205-7533. DOI: 10.1139/f05-179. <URL: http://www.nrcresearchpress.com/doi/abs/10.1139/f05-179> (visited on 07/07/2012).

[2] M. Kéry and M. Schaub. Bayesian population analysis using WinBUGS : a hierarchical perspective. Boston: Academic Press, 2011. ISBN: 9780123870209 0123870208. <URL: http://www.vogelwarte.ch/bpa.html>.

[3] C. J. Krebs. Ecological methodology. 2nd ed. Menlo Park, Calif: Benjamin/Cummings, 1999. ISBN: 0321021738.

[4] M. Plummer. JAGS version 3.3.0 user manual. Oct. 2012. <URL: http://sourceforge.net/projects/mcmc-jags/files/Manuals/3.x/>.

[5] R. C. Team. R: a language and environment for statistical computing. Vienna, Austria, 2013. <URL: http://www.R-project.org>.

Acknowledgements

The organisations and individuals whose contributions have made this analysis report possible include:


blog comments powered by Disqus