Pallid Sturgeon Mortality

Set the number of simulation replicates in years

nReps <- 10000

Time

The variable “doy” is day of year, 1 through 365

doy <- rep(seq(1, 365, by = 1), nReps)
length(doy)

[1] 3650000

head(doy)

[1] 1 2 3 4 5 6

summary(doy)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1      92     183     183     274     365 

The variable “year” is the year of the simulation, starting with 1 and proceeding to nReps.

year <- rep(1:nReps, each = 365)

River stage

The variable “stage” is the river stage in feet. This is simulated using the arima.sim function. The parameters were estimated in a different file [insert file name].

sqrt(2563)

[1] 50.62608

mean(stage)

[1] 383.2532

sd(stage)

[1] 261.314

sd(stage) / mean(stage)

[1] 0.6818313

summary(stage)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -748.6   204.7   383.2   383.3   561.5  1601.0 

Here are the first few years of simulated stage:

plot(stage[1:(365*5)], type = "l", xlab = "Day", ylab = 'Stage (ft)')

Diversion flow rate

The variable “dcfd” is the flow of water through the diversion in cubic feet per day. This is calculated as a linear function of river stage with a minimum and maximum, which assumes a rectangular notch in the levy with passive flow.

This is a conversion factor for converting cubic feet per second to cubic feet per day.

secToDay <- 60*60*24

Since we do not have operational parameters for the diversion, here we are just setting the minimum and maximum flows at the 25th and 75th percentiles of simulated river stage, respectively.

minStage <- quantile(stage, probs = 0.25, na.rm = TRUE)
maxStage <- quantile(stage, probs = 0.75, na.rm = TRUE)
minDivCfs <- 0
maxDivCfs <- 75000
minDivCfd <- minDivCfs * secToDay
maxDivCfd <- maxDivCfs * secToDay

The slope of the relationship between diversion flow rate and river stage is currently arbitrary.

cfsSlope <- (maxDivCfd - minDivCfd) / (maxStage - minStage) # rise over run

In between the minimum and maximum diversion flows, the diversion flow rate is a linear function of river stage.

dcfd <- stage * cfsSlope - (minStage * cfsSlope)

Below the minStage, the diversion flow rate is zero cfs. Above the maxStage, the diversion flowRate is 75,000 cfs.

dcfd[stage < minStage] <- minDivCfd
dcfd[stage > maxStage] <- maxDivCfd

Plot the relationship between diversion flow rate and river stage.

plot(stage[1:1000], dcfd[1:1000], xlab = 'Stage (ft)', ylab = 'Flow rate (Cubic feet per day)')

Plot discharge over time for a subset of data.

plot(dcfd[1:(365*10)], type = "l", xlab = "Day", ylab = 'Flow rate (Cubic feet per day)')

Water temperature

Low average daily water temperature may result in torpor so that finish are not moving and are, consequently, unavailable to be captured by the diversion. The variable “waterTemp” is a binary variable for which a value of “0” indicates that the average daily water temperature is is sufficiently low to induce torpor, and a value of “1” indicates that it is not low enough to induce torpor. This is a function of the day of year based on historical average daily water temperatures. The parameters were estimated in a different file [insert file name].

The variable “firstDay” is the first cumulative day of the year that is on average above the minimum water temperature threshold, and the variable “lastDay” is the last cumulative day of the year that is on average above the minimum water temperature threshold.

firstDay <- 75
lastDay <- 359

Fill the variable with zeros, then replace with 1 for day of year with temperature above the threshold.

waterTemp <- rep(0, length(doy))
waterTemp[doy >= firstDay & doy <= lastDay] <- 1

Mortality

The variable “mortalityRate” is the number of fish captured by the diversion per cubic feet of water diverted.

mortalityRate <- 1 / (9.472*(10^9)) # 1 fish per that many cubic feet per day of diverted water

The variable “m1” is mortality computed for all days of the year, and the variable “m2” is mortality computed with zero mortality for days of the year that are below the torpidity temperature threshold.

m1 <- mortalityRate * dcfd
m2 <- mortalityRate * dcfd * waterTemp

Plot a few years of daily mortality.

plot(1:(365*5), m1[1:(365*5)], type = "l")

plot(1:(365*5), m2[1:(365*5)], type = "l")

mean(m1)

[1] 0.3421026

sd(m1)

[1] 0.2773446

sd(m1) / mean(m1)

[1] 0.8107059

mean(m2)

[1] 0.2819343

sd(m2)

[1] 0.2877488

sd(m2) / mean(m2)

[1] 1.020624

hist(m1, main = 'Ignoring torpor', xlab = 'Fish taken per day', col = 'grey')

hist(m2, main = 'Including torpor', xlab = 'Fish taken per day', col = 'grey')

Results

First we need to aggregate the data annually, getting the total fish taken per year.

m1Annual <- tapply(m1, INDEX = list(year), FUN = sum)
m2Annual <- tapply(m2, INDEX = list(year), FUN = sum)

hist(m1Annual, main = 'Ignoring torpor', xlab = 'Fish taken per year', col = 'grey')

hist(m2Annual, main = 'Including torpor', xlab = 'Fish taken per year', col = 'grey')

These are the results to be used in RAMAS. These are the average number of fish taken by the diversion per year, and the standard deviation in the number of fish taken per year. …for the m1 scenario

mean(m1Annual)

[1] 125.0085

sd(m1Annual)

[1] 44.79841

…and for the m2 scenario

mean(m2Annual)

[1] 102.9118

sd(m2Annual)

[1] 38.02213