--- title: "Centennial Valley Arctic Grayling Adaptive Management Project 2017 Spring Update" date: '`r format(Sys.time(), "%d %B, %Y")`' output: pdf_document --- ```{r fish_data_import, echo=FALSE, warning=FALSE, results='hide', error=FALSE, message=FALSE} # Set directory and load packages ----------------------------------------- rm(list=ls(all=TRUE)) ######################################################################################################## ## Set report year for subsetting data below (i.e., spawning habitat) ## ## Spawning habitat surveys occur every other year; first complete survey took two years (2014, 2015) ## ######################################################################################################## # report.year <- c(2014,2015) report.year <- 2017 # Set working directory to the AnnualReport folder setwd(paste("C:\\Users\\jewarren\\Documents\\Biology\\Fish & Wildlife\\Fish\\Grayling\\Adaptive Management Project\\AnnualReport\\Report Archive\\", report.year, sep="")) ## Load necessary packages library(knitr) library(RMark) library(lme4) library(xtable) library(FSA) library(plotrix) # for histStack library(spatstat) ###**********************************************************************************### ### Arctic grayling and habitat data for fitting competing models of system dynamics ### ###**********************************************************************************### ## Data are stored in a database co-maintained by J.Warren, USFWS, Lakeview, Montana - 406.276.3536 ## and M. Jaeger, MFWP, Dillon, Montana - 406.683.9310. The data frame has a column for each ## variable found in models used to predict Arctic grayling abundance in Red Rock Creek based ## on spawning population estimates and habitat conditions. # Import fish data -------------------------------------------------------- fish.all.years <- read.csv(paste("AG_models_data_",report.year,".csv", sep=""), header=T) #All fish weir count data, 1962 to current. ## grayling data are corrected for detection probability grayling <- fish.all.years[-(1:7),] #1994 to current to include detection corrected grayling numbers ## Create a time variable for both data sets fish.all.years$time <- fish.all.years$year-1961 grayling$time <- grayling$year-1993 ## Make sure values are numeric (necessary due to NAs in the data set) fish.all.years$year <- as.numeric(fish.all.years$year) fish.all.years$grayling <- as.numeric(fish.all.years$grayling) fish.all.years$cutt <- as.numeric(fish.all.years$cutt) grayling$year <- as.numeric(grayling$year) grayling$grayling <- as.numeric(grayling$grayling) ## Use time-series mean for missing values of annual survival ## These will update annually and the influence of the missing values will decline with time grayling$phi[is.na(grayling$phi)] <- mean(na.omit(grayling$phi)) # Import habitat data ------------------------------------------------------ ## Import winter habtitat data and create a point pattern dataset for inverse-distance weighted smoothing of ## dissolved oxygen and water depth. wnt.hab <- read.csv(paste("Winter_habitat_data_",report.year, ".csv", sep=""), header=T) names(wnt.hab) <- c("site","date","wyear","do.0m","do.1m","do.b","temp.0m","temp.1m","temp.b","snow.depth","ice.depth","h2o.depth","x","y") wnt.hab <- wnt.hab[!(wnt.hab$wyear==1997),] #drops surveys from 1997; uncertain quality of data wnt.hab <- wnt.hab[!(is.na(wnt.hab$do.0m)),]; wnt.hab <- wnt.hab[!(is.na(wnt.hab$h2o.depth)),] #drops records with NA values for 0m dissolved oxygen ## Import spawning habitat data (currently includes Red Rock and Elk Springs creeks) spawn.hab <- read.csv(paste("Spawning_habitat_data_", report.year, ".csv", sep=""), header=T) spawn.hab$date <- as.Date(spawn.hab$date, format="%m/%d/%Y") spawn.hab$year <- format(spawn.hab$date, "%Y") spawn.hab$cdate <- format(spawn.hab$date, "%j") spawn.hab$s.fines <- 0 spawn.hab$s.cobble <- 0 ## Subset spawning habitat data to match report year THIS DIFFERS FROM THE ANNUAL REPORT DUE TO SPAWNING HABITAT BEING ## QUANTIFIED AFTER THE SPRING UPDATE IS COMPLETED spawn.hab <- if(report.year %in% c(2013:2017)) { spawn.hab <- subset(spawn.hab, spawn.hab$year %in% c(2014:2015)) } else if (report.year > 2017) {subset(spawn.hab, spawn.hab$year==report.year-1 | spawn.hab$year==report.year-2)} ## Dam data with predicted passage and calculated back-flooded area -------------------------------- dam <- read.csv(paste("Beaver_dam_data_", report.year, ".csv", sep=""), header=T) ## Import reach data for calculating area of suitable habitat below reach.data <- read.csv("Reach_data.csv", header=T) ``` ```{r demographic_rates, echo=F, warning=FALSE} ###**********************************************************************************### ### Demographic rates needed for the grayling balance equation ### ###**********************************************************************************### ### Length-age relationship used to estimate fish age # Length probability key (package FSA)-------------------------------------------------- # Import back-calculated age at length data for developing age at length key # These data are psuedo-replicated and introduce additional error because age of an individual is back-calculated # to prior years based on width of scale annuli. # NEED TO HAVE QUERY INCLUDE YEAR OF CAPTURE ONLY FOR ALL FISH. DO NOT HAVE YEAR FOR BACK-CALCULATED AGES; NAs IN YEAR # COLUMN USED TO SUBSET INDIVIDUALS TO ESTIMATE ANNUAL PROPORTION 3YOs (IS THIS STILL CORRECT? NO 'NA's' FOR YEAR IN .CSV) # INITIAL KEY BUILT USING BACK-CALCULATED AGES FROM 1994 AND 2015 TRAP DATA. ADDED OTHER YEARS OF TRAP DATA TO ESTIMATE # AGES FOR ESTIMATING PROPORTION 3 YEAR OLDS. NEED TO ADD NEW LENGTH DATA EACH YEAR AND COMBINE TRAP AND ELECTROFISHING # DATA ONCE DATABASE IS SET UP. CURRENTLY OVERESTIMATES THE NUMBER OF FISH HANDLED? 1994 CERTAINLY. 2016 1 FISH OFF? FUNCTION # OF USING BACK-CALCULATED VALUES? Have not figured out how to create a dataframe that has 'NA' where an age and length were # back-calculated; need to figure this out for the purposes of not inflating the number of 3YOs in the spawning run. ## Can/should update known-age fish as scale reading occurs bc_key <- read.csv("BC_length_at_age.csv", header=TRUE) key <- read.csv(paste("Age_at_length_key_data_", report.year, ".csv", sep=""), header=TRUE) key$age <- ifelse(key$tl<255, 1, key$age) #fish less than the smallest length category are manually set to 1 year olds; no data in key to do that key <- key[!is.na(key$tl),] key.age <- Subset(key,!is.na(age)) # get fish without NAs in age variable key.length <- Subset(key,is.na(age)) # get fish with NAs in age variable ## Bin individual length records into 10mm length bins for calculating probability key ## 'startcat' is the smallest length category used to build the key based on; 'w' is bin width # ag.lngth.cat <- lencat(~tl,data=key.age,startcat=75,w=10) ag.lngth.cat <- lencat(~tl,data=bc_key, startcat=75,w=10) #added 6/19/2018; includes back-calculated ages 1994, 2013, 2015. Also includes age-at-capture data 2016-2017. Psuedoreplication? ## Frequency of ages by length bin category ag.raw <- with(ag.lngth.cat,table(LCat,age)) ## Calculates probability age key for age by length bin category ag.key <- prop.table(ag.raw,margin=1) key.length1 <- alkIndivAge(ag.key,age~tl,data=key.length) key.all <- rbind(key.age[!is.na(key.age$year),],key.length1) ## Calculate annual estimate of proportion 3YO based on age-length key predicted ages age.lngth.summary <- tapply(key.all$age, list(key.all$year, key.all$age), length) p3YO.df <- data.frame(as.numeric(row.names(age.lngth.summary)),transform(age.lngth.summary[,3]/rowSums(age.lngth.summary, na.rm=T))) names(p3YO.df) <- c("year","p3YO") grayling <- merge(grayling, p3YO.df, all.x=T, by.x="year", by.y="year") # hist(~age,data=key.all,breaks=0:8,xlim=c(1,8),xlab="Age (yrs)",main="",ylim=c(0,2500)) ## Grayling age class size distributions from Mogen 1996 x <- seq(275, 450, length=700) age2 <- dnorm(x,mean=322,sd=10.6) age3 <- dnorm(x,mean=365,sd=10.5) age4 <- dnorm(x,mean=386,sd=6.0) age5 <- dnorm(x,mean=404,sd=7.9) # plot(x,age2,type="l",ylim=c(0,0.08), lwd=2,col="red") # lines(x,age3,lwd=2,col="blue") # lines(x,age4,lwd=2,col="green") # lines(x,age5,lwd=2,col="black") ## The estimated number of recruits for a given year and females t - 3 are calculated by multiplying ## the proportion of a group (female or recruit) 1) estimated using age-length key (recruits) or proportion of ## females observed during electrofishing on the Red Rock Creek trend section by the estimated abundance corrected ## for detection probability. ## Estimated number of grayling recruits (est3YO) each year grayling$est3YO <- grayling$grayling*grayling$p3YO ## Number of female grayling t-3 grayling$ft3 <- NA for (i in 4:length(grayling$time)) grayling$ft3[i] <- na.pass(grayling$grayling[i-3]*0.5) ## Mean female length t-3 grayling$mn.flt3 <- NA for (i in 4:length(grayling$time)) grayling$mn.flt3[i] <- na.pass(grayling$mn.fl[i-3]) ## Use time-series mean for missing values of mean female length ## These will update annually and the influence of the missing values will decline with time grayling$mn.flt3[is.na(grayling$mn.flt3)] <- mean(na.omit(grayling$mn.flt3)) ###***********************************************************************************************### ## Length-fecundity relationship using data from Lund (1974) and Bishop (1971). Lund provided mean ## number of eggs and lengths by length category; Bishop provided data from measured individuals. ## First 3 observations are Lund's. One of Bishop's observations was excluded as an outlier (obs. #13). ###***********************************************************************************************### # Length-fecundity relationship ------------------------------------------- lngth <- c(33.3,37.6,42.7,41.9,38.3,40.1,39.2,43.2,40.1,40.3,36.8,41.3,25.4,38.5,39.1,37.4,40.9,38.0) eggs <- c(5970,9162,9454,9919,12432,10010,8439,15905,9667,10780,7830,11247,7649,7293,6120,8496,12311,6952) alpha.m <-lm(eggs[-13]~lngth[-13]) #relationship between grayling length and fecundity alpha <- function(x) alpha.m$coefficients[1]+alpha.m$coefficients[2]*x # Egg hatchability rate, Beta, taken from Lund (1974) during 1972 and 1973, Elk Lake, MT. B <- c(0.04, 0.12) ``` ```{r spawn_hab_suitability, echo=FALSE, Warnings=FALSE, eval=T} ###**********************************************************************************************### ### Calculate spawning habitat suitability to fit spawning habitat models ### ###**********************************************************************************************### # Spawning habitat suitability -------------------------------------------- ## Drop site 21B; duplicate survey of 21A NOT NECESSARY? 8/21/2018 # spawn.hab <- droplevels(spawn.hab[spawn.hab$site!="21B",]) ## Calculate riffle area; Two sections have continous riffles - Corral Creek to Antelope Creek and Forks (elk springs creek) spawn.hab$area <- spawn.hab$riffle_lngth*spawn.hab$bf_width ## Calculate percent fines (<=2.8) and percent cobble (8-180); row numbers work with first total column = 9, > 180 = column 26 spawn.hab$p.fines <- rowSums(spawn.hab[9:10])/rowSums(spawn.hab[9:23], na.rm=T) spawn.hab$p.cobble <- rowSums(spawn.hab[13:22])/rowSums(spawn.hab[9:23], na.rm=T) ## Function to predict suitability as a function of percent fines ## See Bolker (2008) pg. 90 for explanation of piecewise linear function spawn.hab$s.fines <- ifelse(spawn.hab$p.fines <= 0.1, 1, ifelse(spawn.hab$p.fines<=0.5, 1 - 2.5*(spawn.hab$p.fines - 0.1),0)) spawn.hab$s.cobble <- ifelse(spawn.hab$p.cobble < 0.2, 5*spawn.hab$p.cobble, 1) ## Calculate suitable area by riffle using the two criteria (m^2) spawn.hab$Ats.riffle <- ifelse(spawn.hab$s.fines + spawn.hab$s.cobble < 2, 0, spawn.hab$area) #absolute suitability spawn.hab$Aws.riffle <- spawn.hab$s.fines*spawn.hab$s.cobble*spawn.hab$area #weighted suitability ## Estimate suitable area by site (m^2 per m); DATAFRAME SHOULD HAVE ONLY CURRENT YEAR DATA (update in database query) ####DELETED 'UNLIST' AT FRONT OF DENOMINATOR 3/7/2018 - DON'T THINK IT WAS DOING ANYTHING#### Ats.site <- round(tapply(spawn.hab$Ats.riffle, spawn.hab$site, sum)/tapply(spawn.hab$site_lngth, spawn.hab$site, unique),4) Aws.site <- round(tapply(spawn.hab$Aws.riffle, spawn.hab$site, sum)/tapply(spawn.hab$site_lngth, spawn.hab$site, unique),4) ## Estimate mean suitability at the reach length ###NEED TO ADD REACH AND CREEK TO SITE? THE BELOW CODE DOESN'T WORK; ADDED CREEK TO SPAWN.HAB TO FIX, BUT DIDN'T FIX THE PROBLEM. NOT SURE THE EASIEST WAY TO GET ATS AND AWS BY STREAM 10/18/2017 site <- row.names(Ats.site) # spawn.hab.reach <- data.frame(data.frame(Ats.site[,max(col(Ats.site))], Aws.site[,max(col(Aws.site))], site, stringsAsFactors = F, row.names=NULL), stringsAsFactors = F, row.names=NULL) #This was for when Ats.site and Aws.site had multiple data from multiple years # spawn.hab.reach <- data.frame(data.frame(site=as.factor(site), Ats.site, Aws.site, stringsAsFactors = F, row.names=NULL), stringsAsFactors = F, row.names=NULL) #NOT CERTAIN WHY THERE WERE NESTED DATA.FRAME STATEMENTS; CHNAGED TO SINGLE DATA.FRAME STATEMENT 5/30/2018 spawn.hab.reach <- data.frame(site=as.factor(site), Ats.site, Aws.site, stringsAsFactors = F, row.names=NULL) spawn.hab.reach <- merge(unique(spawn.hab[,c(1:3)]), spawn.hab.reach, by="site", all.x=T, all.y=F) spawn.hab.reach <- spawn.hab.reach[,c(2:3,1,4:length(spawn.hab.reach[1,]))] ## Subset sites from Red Rock, Elk, and Picnic Creeks spawn.hab.rrc <- droplevels(subset(spawn.hab.reach, spawn.hab.reach$creek %in% c("Red Rock Creek", "Elk Creek"))) spawn.hab.rrc <- merge(spawn.hab.rrc, reach.data[,2:3], by="reach", all.x=T) ## Calculate mean suitable habitat by reach, multiply mean by reach length to get estimated total spawning habitat by reach ## Multiply mean reach estimate by 0.1 to convert to hectares (same as mulitplying by 1000 to convert m^2/m to m^2/km and 0.0001 to convert m^2/km to ha) Ats.rrc <- data.frame(row.names(tapply(spawn.hab.rrc$Ats.site, spawn.hab.rrc$reach, mean)), tapply(spawn.hab.rrc$Ats.site, spawn.hab.rrc$reach, mean), row.names = NULL) #mean Ats by reach in Red Rock Creek names(Ats.rrc) <- c("reach","mn.Ats") Ats.rrc <- merge(Ats.rrc, reach.data[,1:3], by="reach", all.x=T) #add reach and reach data to reach Ats means Ats.rrc$total <- round(Ats.rrc$mn.Ats*Ats.rrc$length*0.1,3) #multiply reach mean Ats by reach length for total area of spawning habitat by reach dam.bw <- round(tapply(dam[dam$year==report.year,]$backwaterlngth, dam[dam$year==report.year,]$reach, sum)/1000,3) dam.bw <- data.frame(row.names(dam.bw), dam.bw, row.names=NULL) names(dam.bw) <- c("reach","length.bw") Ats.rrc <- merge(Ats.rrc, dam.bw, all.x=T, by="reach") Ats.rrc$length.bw[is.na(Ats.rrc$length.bw)] <- 0 Ats.rrc$total.corrected <- round(Ats.rrc$mn.Ats*(Ats.rrc$length-Ats.rrc$length.bw)*0.1,3) #subtract back-watered areas from total reach length Ats.rrc$total.corrected <- ifelse(Ats.rrc$total.corrected <0, 0, Ats.rrc$total.corrected) Aws.rrc <- data.frame(row.names(tapply(spawn.hab.rrc$Aws.site, spawn.hab.rrc$reach, mean)), tapply(spawn.hab.rrc$Aws.site, spawn.hab.rrc$reach, mean), row.names = NULL) #mean Aws by reach in Red Rock Creek names(Aws.rrc) <- c("reach","mn.Aws") Aws.rrc <- merge(Aws.rrc, reach.data[,1:3], by="reach", all.x=T) #add reach and reach data to reach Aws means Aws.rrc$total <- round(Aws.rrc$mn.Aws*Aws.rrc$length*0.1,3) #multiply reach mean Aws by reach length for total area of spawning habitat by reach Aws.rrc <- merge(Aws.rrc, dam.bw, all.x=T, by="reach") Aws.rrc$length.bw[is.na(Aws.rrc$length.bw)] <- 0 Aws.rrc$length.bw[is.na(Aws.rrc$length.bw)] <- 0 Aws.rrc$total.corrected <- round(Aws.rrc$mn.Aws*(Aws.rrc$length-Aws.rrc$length.bw)*0.1,3) #subtract back-watered areas from total reach length Aws.rrc$total.corrected <- ifelse(Aws.rrc$total.corrected <0, 0, Aws.rrc$total.corrected) ## Calculate Ats and Aws for each reach SOMEWHAT REDUNDANT WITH CODE ABOVE FOR RRC (BELOW IS FOR WHOLE DATAFRAME) ## Divided suitable riffle area by site length; summed that by site; multiplied by reach length (km) and by 0.1 to convert to ha ## THIS DOESN'T CURRENTLY ACCOUNT FOR YEAR, ALTHOUGH THERE ARE MULTIPLE YEARS OF DATA, NOR DOES IT ADD YEAR TO THE DATA FRAME - STILL NEED TO DO THAT ## CURRENT SOLUTION IS TO USE THE QUERY TO SELECT THE CURRENT REPORT YEAR'S DATA # Ats.reach <- tapply(spawn.hab.reach$Ats.site, spawn.hab.reach$reach, sum)*tapply(spawn.hab.reach$reach.lngth, spawn.hab.reach$reach, mean)*0.1 # Aws.reach <- tapply(spawn.hab$Aws.riffle, spawn.hab$reach, sum)*tapply(spawn.hab$reach.lngth, spawn.hab$reach, mean)*0.1 # write.csv(spawn.hab, "spawn.hab.csv", row.names=F) Riffle_suitability <-data.frame(I(spawn.hab$site), spawn.hab$riffle, spawn.hab$Ats.riffle, spawn.hab$Aws.riffle) names(Riffle_suitability) <- c("site","riffle","Ats.riffle","Aws.riffle") # write.csv(Riffle_suitability, "Riffle_suitability.csv", row.names=F) # write.csv(cbind(Ats.reach,Aws.reach), "Reach_suitability.csv") ## NEED TO ADD ANNUAL ESTIMATES OF Ats and Aws TO GRAYLING DATAFRAME ``` ```{r spawn_hab_suitability_OLD, echo=FALSE, Warnings=FALSE, eval=F} ###**********************************************************************************************### ### Calculate spawning habitat suitability to fit spawning habitat models ### ###**********************************************************************************************### # Spawning habitat suitability -------------------------------------------- ## Drop site 21B; duplicate survey of 21A spawn.hab <- droplevels(spawn.hab[spawn.hab$site!="21B",]) ## Calculate riffle area; Two sections have continous riffles - Corral Creek to Antelope Creek and Forks (elk springs creek) spawn.hab$area <- spawn.hab$riffle_lngth*spawn.hab$bf_width ## Calculate percent fines (<=2.8) and percent cobble (8-180); row numbers work with first total column = 9, > 180 = column 26 spawn.hab$p.fines <- rowSums(spawn.hab[9:10])/rowSums(spawn.hab[9:23], na.rm=T) spawn.hab$p.cobble <- rowSums(spawn.hab[13:22])/rowSums(spawn.hab[9:23], na.rm=T) ## Function to predict suitability as a function of percent fines ## See Bolker (2008) pg. 90 for explanation of piecewise linear function spawn.hab$s.fines <- ifelse(spawn.hab$p.fines<=0.1, 1, ifelse(spawn.hab$p.fines<=0.5, 1 - 2.5*(spawn.hab$p.fines - 0.1),0)) spawn.hab$s.cobble <- ifelse(spawn.hab$p.cobble<0.2, 5*spawn.hab$p.cobble, 1) ## Calculate suitable area by riffle using the two criteria (m^2) spawn.hab$Ats.riffle <- ifelse(spawn.hab$s.fines + spawn.hab$s.cobble < 2, 0, spawn.hab$area) #absolute suitability spawn.hab$Aws.riffle <- spawn.hab$s.fines*spawn.hab$s.cobble*spawn.hab$area #weighted suitability # ## Calculate suitable area by riffle using the two criteria COMMENTED OUT 6/19/18; WASN'T WORKING # spawn.hab$Ats.riffle <- ifelse(spawn.hab$s.fines+spawn.hab$s.cobble < 2, 0, spawn.hab$area/spawn.hab$site.lngth) #absolute suitability # spawn.hab$Aws.riffle <- ((spawn.hab$s.fines*spawn.hab$area + spawn.hab$s.cobble*spawn.hab$area)/2)/spawn.hab$site.lngth #weighted suitability # write.csv(spawn.hab, "spawn.hab.csv", row.names=F) Riffle_suitability <- as.data.frame(cbind(spawn.hab$site, spawn.hab$riffle, spawn.hab$Ats.riffle, spawn.hab$Aws.riffle)) names(Riffle_suitability) <- c("site","riffle","Ats.riffle","Aws.riffle") # write.csv(Riffle_suitability, "Riffle_suitability.csv", row.names=F) ## Calculate Ats and Aws based on area of riffles ## Divided suitable riffle area by site length; summed that by site; multiplied by reach length (km) and by 0.1 to convert to ha Ats.reach <- tapply(spawn.hab$Ats.riffle, spawn.hab$reach, sum)*tapply(spawn.hab$reach_lngth, spawn.hab$reach, mean)*0.1 Aws.reach <- tapply(spawn.hab$Aws.riffle, spawn.hab$reach, sum)*tapply(spawn.hab$reach_lngth, spawn.hab$reach, mean)*0.1 # write.csv(cbind(Ats.reach,Aws.reach), "Reach_suitability.csv") ## NEED TO ADD ANNUAL ESTIMATES OF Ats and Aws TO GRAYLING DATAFRAME ``` ```{r fish_passage, echo=FALSE, warnings=FALSE, eval=FALSE} ###**********************************************************************************************### ### Calculate probability of passage and area flooded for beaver dams ### # NOT CURRENTLY DOING THIS BECAUSE DAMS ARE BEING BREACHED 2017-2020 # Add reach name to dam data.frame NOT NECESSARY AS OF 8/21/2018, REACH ALREADY IN SPREADSHEET # dam$reach <- cut(dam$dist_lake, breaks=c(rev(reach.data$b_dist_lk), max(reach.data$t_dist_lk)), labels=rev(reach.data$reach)) # Add a column of cumulative passage probability to the dam data.frame dam$cpass.prob <- unlist(tapply(dam$p_pass, dam$year, cumprod)) dam$order <- as.numeric(row.names(dam)) # Add reach top and bottom distances dam <- merge(dam, reach.data[,c(2,8,9)], by="reach", sort=F) dam <- dam[order(dam$order),] # Minimum distance to lake for each year min.yr <- data.frame(year=2013:max(dam$year), dist.lake=rep(NA,length(2013:max(dam$year)))) for (i in 1:length(unique(dam$year))) min.yr[i,2] <- min(na.exclude(dam$dist_lake[dam$year==(i + min(unique(dam$year))-1)])) ``` ```{r wint_hab_suitability, echo=FALSE, warnings=FALSE} ###*********************************************************************************************### ### Calculate winter habitat suitability to fit winter habitat models ### ###*********************************************************************************************### # Winter habitat suitability ---------------------------------------------- do.mn <- function(x,y) rowMeans(cbind(x,y), na.rm=TRUE) ## Addressed difference between Gangloff and Davis methods (Gangloff measured top and bottom of h2o column, Davis measured top [0m], 1m, and 2 m if h2o >1.9m) wnt.hab$do <- ifelse(!(is.na(wnt.hab$do.0m)) & !(is.na(wnt.hab$do.1m)), do.mn(wnt.hab$do.0m,wnt.hab$do.1m), ifelse(!(is.na(wnt.hab$do.0m)) & !(is.na(wnt.hab$do.b)) & wnt.hab$h2o.depth>=1 & wnt.hab$h2o.depth<1.5, do.mn(wnt.hab$do.0m,wnt.hab$do.b), ifelse(!(is.na(wnt.hab$do.0m)) & !(is.na(wnt.hab$do.b)) & wnt.hab$h2o.depth>=1.5, do.mn(wnt.hab$do.0m,do.mn(wnt.hab$do.0m,wnt.hab$do.b)), wnt.hab$do.0m))) wnt.hab$date <- as.Date(wnt.hab$date, format="%m/%d/%Y") wnt.hab$year <- as.numeric(format(wnt.hab$date, "%Y")) wnt.hab$month <- as.numeric(format(wnt.hab$date, "%m")) wnt.hab$smp.period <- as.factor(paste(wnt.hab$month, wnt.hab$wyear, sep="-")) wnt.hab <- subset(wnt.hab, wnt.hab$month != 4) #excludes April 2011 data; suspiciously low values and few sites sampled ## Import Upper Red Rock Lake polygon vertices urrl.poly <- read.csv("URRL_vertices.csv", header=T) urrl.poly[1,2] <- 440500 #change the first x value because the polygon can not be closed ## Add columns to match 'wnt.hab.subset' below to allow an rbind for using shore vertices in calculating depth (shore = 0) urrl.poly$site <- NA urrl.poly$date <- NA urrl.poly$wyear <- NA urrl.poly$do.0m <- NA urrl.poly$do.1m <- NA urrl.poly$do.b <- NA urrl.poly$temp.0m <- NA urrl.poly$temp.1m <- NA urrl.poly$temp.b <- NA urrl.poly$snow.depth <- NA urrl.poly$ice.depth <- NA urrl.poly$h2o.depth <- 0 urrl.poly$do <- 0 urrl.poly$year <- NA urrl.poly$month <- NA urrl.poly$smp.period <- NA urrl.poly <- urrl.poly[,c(4:15,2,3,16:19)] ## Use sampling period to create a function to calculate suitable habitat for each sampling period and adding 0 water depth at Upper Lake ## boundary points. Steps are 1) calculate suitable area by sampling occassion, 2) select minimum area (or average?) for each year, 3) put value ## into working dataframe for fitting models ## Set the window for calculating inverse-distance weighted values of DO and water depth; for a square window xrange=c(440000,445250) yrange=c(4938000, 4942500) wnt.hab.poly.data <- data.frame(x=rev(urrl.poly$x), y=rev(urrl.poly$y)) xrange=c(min(urrl.poly$x),max(urrl.poly$x)) yrange=c(min(urrl.poly$y),max(urrl.poly$y)) wnt.hab.window <- owin(xrange=xrange, yrange=yrange, poly=wnt.hab.poly.data, unitname="meters") ## 1) Subset winter habitat data by sampling period, 2) create a ppp object for each subset, 3) interpolate DO and depth for each subset, and ## 4) extract area of suitable winter habitat for each subset. Results are stored as lists on the workspace wnt.hab.subset <- list() wnt.hab.ppp <- list() wnt.hab.idw <- list() wnt.hab.suitable <- list() wnt.hab.suitable.do <- list() wnt.hab.suitable.depth <- list() wnt.do.area <- vector() wnt.depth.area <- vector() wnt.hab.area <- vector() for (i in 1:length(unique(wnt.hab$smp.period))) { wnt.hab.subset[[i]] <- wnt.hab[wnt.hab$smp.period==as.vector(unique(wnt.hab$smp.period))[i],] wnt.hab.subset[[i]] <- rbind(wnt.hab.subset[[i]], urrl.poly) wnt.hab.subset[[i]]$smp.period <- wnt.hab.subset[[i]]$smp.period wnt.hab.subset[[i]]$wyear <- wnt.hab.subset[[i]]$wyear wnt.hab.ppp[[i]] <- ppp(x=wnt.hab.subset[[i]]$x, y=wnt.hab.subset[[i]]$y, n=1, window=wnt.hab.window, marks=data.frame(wnt.hab.subset[[i]]$do, wnt.hab.subset[[i]]$h2o.depth), unitname="meters") wnt.hab.idw[[i]] <- idw(wnt.hab.ppp[[i]], power=2, at="pixels") wnt.hab.suitable[[i]] <- solutionset(wnt.hab.idw[[i]][[1]] >= 4 & wnt.hab.idw[[i]][[2]] >= 100) wnt.hab.suitable.do[[i]] <- solutionset(wnt.hab.idw[[i]][[1]] >= 4) wnt.hab.suitable.depth[[i]] <- solutionset(wnt.hab.idw[[i]][[2]] >= 100) wnt.do.area[i] <- area(wnt.hab.suitable.do[[i]]) wnt.depth.area[i] <- area(wnt.hab.suitable.depth[[i]]) wnt.hab.area[i] <- area(wnt.hab.suitable[[i]]) } ## Select the minimum annual suitable habitat estimate based on water year smp.period <- as.vector(unique(wnt.hab$smp.period)) wyear <- vector() for (i in 1:length(unique(wnt.hab$smp.period))) wyear[i] <- as.numeric(strsplit(smp.period, "-")[[i]][2]) wnt.hab.area <- data.frame(wyear,smp.period,wnt.hab.area) names(wnt.hab.area)[3] <- "area" wnt.hab.min <- tapply(wnt.hab.area$area, wnt.hab.area$wyear, min) wnt.hab.min <- cbind(as.numeric(row.names(data.frame(wnt.hab.min))), data.frame(wnt.hab.min, row.names=NULL)) names(wnt.hab.min) <- c("year","Wt") #year is water year (Oct. 1 - Sep. 30) grayling <- merge(grayling, wnt.hab.min, all.x=T, by="year") ``` ```{r winter_prediction, echo=FALSE} wnt.hab.pred.fn <- function (N0, Wt, a, b) { Nt <- matrix(NA, nrow=4, ncol=3) #build a 3 x t matrix to hold predicted values (N, P, W) Nt[1:3,1] <- N0 #input initial population values from the grayling data #Nt[,4] <- exp(rnorm(4) * 0.1 - (0.1^2)/2) #create random observation error - NOT CURRENTLY USED # Predict Nt based on the winter habitat model Nt[1:3,3] <- (Wt+5)/Nt[1:3,1] # Estimated available winter habitat Nt[1:3,2]<-a*Nt[1:3,3]/(b + Nt[1:3,3]) # Predict proportional change in survival (Nt[,2]) for each observed value of Wt and set of values for a and b Nt[4,1] <- Nt[3,1]*0.74*Nt[3,2] + Nt[1,1]*0.5*8628.732*0.08*0.014*Nt[1,2]*0.48*Nt[2,2]*0.68*Nt[3,2]*0.87 # Predict grayling abundance based on the winter habitat model return(Nt) } ## Current-year's report prediction for updating model weights; uses the three year's grayling estimates prior to the current year (e.g., 2013, 2014, and 2015 for 2016 prediction) and winter habitat for the three most recent winters (e.g., 2014, 2015, and 2016 for 2016 prediction). This uses the prior year's grayling estimate as the number of fish entering the winter period, and the water year's winter habitat estimate to predict spawning runs. wnt.pred <- wnt.hab.pred.fn(grayling$grayling[(length(grayling$grayling)-3):(length(grayling$grayling)-1)], wnt.hab.min$Wt[(length(wnt.hab.min$Wt)-2):length(wnt.hab.min$Wt)]/10000,1,0.185) ## Last year's prediction wnt.pred.t_1 <- wnt.hab.pred.fn(grayling$grayling[(length(grayling$grayling)-4):(length(grayling$grayling)-2)], wnt.hab.min$Wt[(length(wnt.hab.min$Wt)-2):length(wnt.hab.min$Wt)]/10000,1,0.185) ``` ```{r spawning_hab_prediction, echo=FALSE} spw.hab.pred.fn <- function (N0, Ht, a, b) { Nt <- matrix(NA, nrow=4, ncol=4) #build a 4 x t matrix to hold predicted values (N, P, W) Nt[1:3,1] <- N0 #input initial population values from the grayling data (1994 & 1995) #Nt[,4] <- exp(rnorm(4) * 0.1 - (0.1^2)/2) #create random observation error # Predict Nt based on the spring habitat model Nt[1:3,3] <- Ht/Nt[1:3,1] # Estimated available spawning habitat Nt[1:3,2]<-a*Nt[1:3,3]/(b + Nt[1:3,3]) # Predict recruitment (Nt[,2]) for each observed value of Ht and set of values for a and b Nt[4,1] <- Nt[3,1]*0.52 + Nt[1,1]*0.5*8628.732*Nt[1,2]*0.25*0.44*0.74 # Predict grayling abundance based on the spawning habitat model; multiply adults by *Nt[3,4] for random observer error return(Nt) } ## Current-year's report prediction for updating model weights; uses the three years prior to the current year (e.g., 2013, 2014, and 2015 for 2016 prediction). THIS WILL BE OFF ONE YEAR UNTIL WINTER SAMPLING IS COMPLETED THE NEXT WATER YEAR. SET UP THIS WAY SO WE CAN UPDATE MODELS AND GET THE NEXT YEAR'S PREDICTION IN THE ANNUAL REPORT spw.pred <- spw.hab.pred.fn(grayling$grayling[(length(grayling$grayling)-3):(length(grayling$grayling)-1)],rep(sum(Aws.rrc$total.corrected)/2,3),0.0042,0.0217) ## Next year's prediction spw.pred.t1 <- spw.hab.pred.fn(grayling$grayling[(length(grayling$grayling)-2):length(grayling$grayling)], rep(sum(Aws.rrc$total.corrected)/2,3),0.0042,0.0217) ## Last year's prediction spw.pred.t_1 <- spw.hab.pred.fn(grayling$grayling[(length(grayling$grayling)-4):(length(grayling$grayling)-2)], rep(sum(Aws.rrc$total.corrected)/2,3),0.0042,0.0217) ``` ```{r nonnative_prediction, echo=FALSE} nonnative.pred.fn <- function (N0, Ct, a, b) { Nt <- matrix(NA, nrow=4, ncol=4) #build a 4 x t matrix to hold predicted values (N, P, W) Nt[1:3,1] <- N0 #input initial population values from the grayling data (1994 & 1995) #Nt[,4] <- exp(rnorm(4) * 0.1 - (0.1^2)/2) #create random observation error # Predict Nt based on the spring habitat model Nt[1,3] <- Ct # Estimated number of cutthroat trout in Upper Lake during a cohort's first two years Nt[1:3,2]<-a*Nt[1:3,3]/(b + Nt[1:3,3]) # Predict recruitment (Nt[,2]) for each observed value of Ht and set of values for a and b Nt[4,1] <- Nt[3,1]*0.52 + Nt[1,1]*0.5*8628.732*0.08*(1-Nt[1,2])*0.74 # Predict grayling abundance based on the nonnative predation model #Nt[4,1] <- Nt[3,1]*0.52*Nt[3,4] + Nt[1,1]*0.5*8628.732*0.08*(1-Nt[1,2])*0.74 # Predict grayling abundance based on the nonnative predation model with observer error return(Nt) } ## Current-year's report prediction for updating model weigths; uses the three years prior to the current year (e.g., 2013, 2014, and 2015 grayling estimates and 2014 and 2015 cuttthroat estimates (for preceding winter abundance) for 2016 prediction). nonnative.pred <- nonnative.pred.fn(grayling$grayling[(length(grayling$grayling)-3):(length(grayling$grayling)-1)],mean(grayling$cutt[(length(grayling$grayling)-2):(length(grayling$grayling)-1)]),0.999,1.225) ## Next year's prediction nonnative.pred.t1 <- nonnative.pred.fn(grayling$grayling[(length(grayling$grayling)-2):length(grayling$grayling)], mean(grayling$cutt[(length(grayling$grayling)-1):length(grayling$grayling)]),0.999,1.225) ## Last year's prediction nonnative.pred.t_1 <- nonnative.pred.fn(grayling$grayling[(length(grayling$grayling)-4):(length(grayling$grayling)-2)],mean(grayling$cutt[(length(grayling$grayling)-3):(length(grayling$grayling)-2)]),0.999,1.225) ``` ```{r model_weight_updating, echo=FALSE} ## Update priors each year using the prior year's annual report UNLIKE THE ANNUAL REPORT .Rmd, THIS REPORT DOES NOT UPDATE MODEL PRIORS FOR THE SUBSEQUENT YEAR priors <- read.csv("model.wts.csv") #2016 had uninformative priors - c(0.333,0.333,0.333) wnt.prob <- dnorm(grayling$grayling[grayling$year==max(grayling$year)], mean=round(wnt.pred[4,1],0), sd=round(sd(grayling$grayling, na.rm=T),0)) spw.prob <- dnorm(grayling$grayling[grayling$year==max(grayling$year)], mean=round(spw.pred[4,1],0), sd=round(sd(grayling$grayling, na.rm=T),0)) nonnative.prob <- dnorm(grayling$grayling[grayling$year==max(grayling$year)], mean=round(nonnative.pred[4,1],0), sd=round(sd(grayling$grayling, na.rm=T),0)) model.probs <- c(wnt.prob,spw.prob,nonnative.prob) model.wts <- priors*model.probs/sum(priors[1,1]*model.probs[1],priors[2,1]*model.probs[2],priors[3,1]*model.probs[3]) mdl.wtd.pred <- round(model.wts[1,1]*wnt.pred[4,1]+model.wts[2,1]*spw.pred[4,1]+model.wts[3,1]*nonnative.pred[4,1],0) ``` * The Centennial Valley Arctic Grayling Adaptive Management Plan (AMP) is being implemented to identify limiting factor(s) for Arctic grayling in the upper Centennial Valley (CV) of southwestern Montana. Non-native hybrid Yellowstone cutthroat trout, spawning habitat, and overwinter habitat have been identified as the three most likely factors that could limit long-term viability of grayling in the upper CV. Long-term viability is expected to be maintained by 1) conserving genetic diversity, 2) establishing spawning and/or refugia in at least two tributaries, and 3) maintaining a spawning population of $\geq$ 1,000 fish. The latter is based on the Species Status Assessment Workshop for Arctic grayling conducted in 2014. * The AMP focus is on identifying factors that cause the spawning population to decline below 1,000 fish and, if that occurs, management actions that will most effectively return the population to objective. An emphasis on learning through 'management as experiment' during the first phase of the AMP is being accomplished via two experiments that 1) reduced non-native hybrid Yellowstone cutthroat trout population (2013--2016) and 2) will maximize availability of spawning habitat (2017--2020). Experiments are being developed for altering winter habitat. To date, natural variability has provided opportunity to explore the hypothesized relationship between grayling spawning population and area of suitable winter habitat in Upper Red Rock Lake (Upper Lake). * The estimated number of Arctic grayling in the `r max(grayling$year)` Red Rock Creek spawning population was **`r grayling$grayling[grayling$year==max(grayling$year)]` (95% CI = `r grayling$lci[grayling$year==max(grayling$year)]`--`r grayling$uci[grayling$year==max(grayling$year)]`**; Figure 1). * The estimated number of Yellowstone cutthroat trout in the Red Rock Creek spawning population was **`r grayling$cutt[grayling$year==max(grayling$year)]`**, an approximate reduction of **`r round(1-grayling$cutt[grayling$year==max(grayling$year)]/grayling$cutt[grayling$year==2014],2)*100`%** from the highest estimated population in 2014 ($\hat{N}$ = `r grayling$cutt[grayling$year==2014]`; Figure 1). Estimated angler harvest was 176 trout. ```{r grayling_cutt_fig, echo=FALSE, fig.align='center', fig.width=6, fig.height=3} ## Plot grayling and cutthroat estimates 2013-current par(mfrow=c(1,2), oma=c(0,0,0,0), mar=c(4,4,3,1), mgp=c(2,1,0)) #mar=c(4,4,2,1), mgp=c(2,1,10)) plot(grayling$grayling[20:length(grayling$grayling)]~grayling$year[20:length(grayling$grayling)], pch=16, ylab="Abundance", xlab="Year", main=paste("a) ", "2013-", max(grayling$year), sep=""), xaxp=c(2013, max(grayling$year), max(grayling$year)-2013), ylim=c(0,3500), cex=.75, cex.axis=.75, cex.lab=.75, cex.main=.75) segments(grayling$year, grayling$lci, grayling$year, grayling$uci) points(grayling$cutt[20:length(grayling$grayling)]~grayling$year[20:length(grayling$grayling)], pch=8, cex=0.75) legend('topright', c("Arctic grayling", "Cutthroat trout"), cex=0.75, pch=c(16,8), bty="n") hist(key.all$age[key.all$year==max(key.all$year)], breaks=c(seq(0.5,8.5,1)), main=paste("b)", max(key.all$year),"Age Distribution"), ylab="Frequency", xlab="Age", cex.axis=.75, cex.lab=.75, cex.main=.75) ``` **Figure 1.** a) Arctic grayling and non-native hybrid Yellowstone cutthroat trout abundance estimates and 95% confidence intervals (grayling only) from Red Rock Creek, 2013--`r max(grayling$year)`, and b) age distribution of the `r max(grayling$year)` grayling spawning population. \pagebreak * Two grayling ($\approx$ 1% of the spawning population) were estimated as caught and released by anglers during angler access surveys conducted mid-April--mid-May; no grayling were reported as caught based on angler self-reporting catch cards. Fourteen of 698 grayling (2%) were observed with hook scars at the fish weir in 2015, the only year these data are available. * We will continue to learn how grayling respond to Yellowstone cutthroat trout population reductions, the first management experiment undertaken as part of the AMP, as 1) grayling cohorts spawned during low trout abundances recruit, and 2) Yellowstone cutthroat trout spawning population recovers. * Suitable winter habitat within Upper Lake (i.e., water depth below the ice $\geq$ 1 m and dissolved oxygen $\geq$ 4 ppm) reached a minimum during February sampling at an estimated `r round(grayling$Wt[grayling$year==max(grayling$year)]/10000, 0)` ha. Grayling spawning population was reduced to $\leq$ 214 fish in all years when <10 ha of suitable winter habitat was available in Upper Lake (Figure 2). * Suitable spawning habitat was most recently quantified in `r max(spawn.hab$year)`, with an estimated total area of suitable spawning habitat ($A_{ts}$) of `r round(sum(Ats.rrc$total.corrected),1)` ha, and weighted area of suitable habitat ($A_{tw}$) of `r round(sum(Aws.rrc$total.corrected),1)` ha, in Red Rock and Elk Springs creeks. Surveys to estimate area of suitable spawning habitat will be completed again this year. ```{r winter_hab_fig, echo=FALSE, fig.height=3, fig.width=7.5, fig.align='center'} ## Example plot of predicted area of suitable winter habitat in Upper Lake for the prior winter - NEEDS ANNUAL UPDATING ONCE MINIMUM HABITAT IS DETERMINED (FOUND IN wnt.hab.area; USE ROW NUMBER FOR LOWEST ESTIMATE FOR THE CURRENT UPDATE YEAR) par(mfrow=c(1,3), oma=c(0,0,0,0), mar=c(4,4,3,1), mgp=c(2,1,0), bty="n") #=c(1,1,2,1), bty="n") plot(wnt.hab.suitable[[19]], main="a", xaxt=FALSE, bty="n") #UPDATE ANNUALLY TO CORRECT ROW WITH MINIMUM ANNUAL WINTER HABITAT ESTIMATE lines(wnt.hab.poly.data, lwd=2) arrows(444500, 4941000, 444500, 4941800, code=2, lwd=4) text(444500, 4942150, "N", font=2, cex=1.5) text(443250, 4939000, "Upper Red Rock Lake", font=2) text(443250, 4938700, paste(round(min(wnt.hab.area$area[wnt.hab.area$wyear==max(wnt.hab.area$wyear)])/10000, 0), "ha", sep=" "), font=2) ## Plot of Minimum Suitable Winter Habitat (ha) by water year plot(wnt.hab.min$Wt/10000, pch=16, main="b", ylab="Suitable Winter Habitat (ha)", xlab="Water Year", ylim=c(0,300), xaxt="n") text(wnt.hab.min$Wt/10000, labels=c(wnt.hab.min$year), pos=c(4, rep(3, length(wnt.hab.min$year)-2), 2), offset=0.25, cex=0.75) # axis(1, labels=c(wnt.hab.min$year), at=1:length(wnt.hab.min$year)) axis(1, labels=NA, at=1:length(wnt.hab.min$year)) plot(grayling$grayling~round(grayling$Wt/10000,0), pch=16, main="c", ylab="Arctic Grayling Spawning Population", xlab="Suitable Winter Habitat (ha)") text(x=grayling[grayling$Wt>180000 & !(is.na(grayling$grayling)),]$Wt/10000, y=grayling[grayling$Wt>180000 & !(is.na(grayling$grayling)),]$grayling, labels=c(grayling[grayling$Wt>180000 & !(is.na(grayling$grayling)),]$year), pos=c(4,4,4), offset=0.25, cex=0.75) rect(xleft=10, ybottom=0, xright=25, ytop=2500, border=NA, col='lightgray') #text(grayling$Wt[grayling$year %in% c(1995,2011,2014:max(grayling$year))]/10000, grayling$grayling[grayling$year %in% c(1995,2011,2014:max(grayling$year))], labels=row.names(sort(wnt.hab.min$Wt))[-c(2,8)], pos=c(4, rep(3, length(wnt.hab.min$year)-2), 2), offset=0.25, cex=0.6) points(grayling[grayling$year %in% c(2016,2017,2018),]$grayling~round(grayling[grayling$year %in% c(2016,2017,2018),]$Wt/10000,0), pch=16) #put points in lower left quadrant back on top of gray rectangle ``` **Figure 2.** a) Extent of minimum area of suitable Arctic grayling winter habitat in Upper Red Rock Lake, `r max(grayling$year)`, b) annual estimate of minimum area of suitable habitat for water years `r min(wnt.hab.area$wyear)`--`r max(wnt.hab.area$wyear)`, and c) grayling spawning population as a function of minimum area of suitable winter habitat for years when both were estimated (1995 [0 ha], 2016 [8 ha], and 2017 [9 ha] points are plotted but not labelled). The shaded polygon represents an hypothesized threshold (10--25 ha) of suitable winter habitat where 1) enough winter habitat is available to sustain grayling population at objective (*N* $\geq$ 1,000 fish, > 25 ha suitable habitat), and 2) winter habitat presumably reduces grayling survival, resulting in grayling population below objective ($\hat{N} \leq$ 214, <10 ha suitable habitat). \pagebreak * The *Winter Habitat*, *Spawning Habitat*, and *Non-native Fish* models predicted `r round(wnt.pred[4,1],0)`, `r round(spw.pred[4,1],0)`, and `r round(nonnative.pred[4,1],0)` grayling, respectively, in the `r max(grayling$year)` Red Rock Creek spawning population. The *Winter Model* continues to be the most supported model, although the *Spawning Habitat Model* was also relatively well supported (Table 1). The *Non-native Fish Model* continued to poorly predict grayling population. Thus, data collected to date provide support for overwinter habitat in Upper Lake and spawning habitat being drivers of grayling abundance in the upper CV. ```{r model_table, echo=F} model_table <- data.frame("Model" = c("Winter Habitat","Spawning Habitat", "Non-native Fish"), "CYPrediction" = c(round(wnt.pred[4,1],0), round(spw.pred[4,1],0), round(nonnative.pred[4,1],0)), "Observed" = c(rep(grayling$grayling[grayling$year==max(grayling$year)], 3)), "Weights" = round(model.wts, 3)) ``` **Table 1.** Arctic grayling spawning abundance model predictions, observed abundance, and relative model weights for `r max(grayling$year)`. Model weights, which sum to 1, are a measure of relative support for a model given the data. ```{r table_2, echo=F, results='asis'} table2 <- xtable(model_table, align="llccc", digits=c(0,0,0,0,3)) names(table2)[2:4] <- c(paste(max(grayling$year), "Prediction"), "Observed", "Model Weights") print(table2, table.placement="H", comment=F, include.rownames=F, floating=FALSE) ``` * Management actions to assess the effect of spawning habitat on grayling abundance were intiated this spring. Spawning habitat was maximized by 1) restoring connectivity to, and habitat within, Elk Springs Creek, and 2) providing access to all spawning habitat in Red Rock Creek by breaching beaver dams (*n* = 14). However, the influence of Yellowstone cutthroat trout abundance and spawning habitat on grayling recovery from the recent decline is currently confounded because the system has both low abundance of Yellowstone cutthroat trout and high per-capita suitable spawning habitat. Therefore, strong grayling cohorts produced by the 2017 or 2018 spawning populations, as predicted by both the *Spawning Habitat* and *Non-native Fish* models, could be related to either, or both, low Yellowstone cutthroat trout abundance or high per-capita spawning habitat. Disentangling the relative influence of these two factors requires maintaining one relatively constant while allowing the other to vary. This is currently being achieved by the second management experiment of the AMP, which maximizes the per-capita area of spawning habitat available to grayling until Yellowstone cutthroat trout return to relatively high abundances (i.e., as soon as 2020). * No management actions to improve winter habitat are presently identified or planned. However, the recent restoration of Elk Springs Creek to improve spawning habitat may also improve winter habitat in Upper Red Rock Lake by allowing the creek to largely circumvent Swan Lake, increasing the flow of highly-oxygenated creek water into Upper Lake. Assessing the potential for Elk Springs Creek restoration, or future management actions, to mitigate winter habitat is currently hampered by low grayling abundance. Based on three years where winter habitat and grayling spawning population were estimated, a threshold level of 10--25 ha of winter habitat appears necessary to overwinter grayling populations greater than the 1,000 fish objective (Figure 2). Based on winter habitat surveys conducted in 2017, the restoration may not have increased area of suitable habitat as currently defined, although during periods of low levels of dissolved oxygen in Upper Lake it may still provide refugia with suitable dissolved oxygen in water depths less than currently considered suitable by the *Winter Habitat* model. ##LITERATURE CITED Boyd, J. 2008. Effects of water temperature and angling on mortality of salmonids in Montana streams. Thesis. Montana State University, Bozeman, Montana. 56 pgs. \pagebreak ###Authors: Jeff Warren$^{1}$, Matt Jaeger$^{2}$, Andrew Gilham$^{3}$, Kyle Cutting$^{4}$, Lucas Bateman$^{2}$, Terrill Paterson$^{2}$. $^{1}$ U.S. Fish and Wildlife Service, Division of Science Resources, $^{2}$ Montana Fish, Wildlife and Parks, $^{3}$ U.S. Fish and Wildlife Service, Montana Fish and Wildlife Conservation Office, $^{4}$ U.S. Fish and Wildlife Service, Red Rock Lakes National Wildlife Refuge.