-
Poulet Camille authored30f3ba4a
#----------------------------------------------
# Temperature effect
#----------------------------------------------
# temperature effect in GR3D
# from Rosso
temperatureEffect = function(tempWater, Tmin, Topt, Tmax){
response = (tempWater - Tmin) * (tempWater - Tmax) / ((tempWater - Tmin) * (tempWater - Tmax) - (tempWater - Topt)^2)
response[tempWater <= Tmin | tempWater >= Tmax] = 0
return(response)
}
thermalRange = function(pct = 0.8, Tmin, Topt, Tmax){
lower = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax) - pct,
interval = c(Tmin, Topt))$root
upper = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax) - pct,
interval = c(Topt, Tmax))$root
return(c(lower = lower,upper = upper))
}
#optimalThermalRange(Tmin = 3, Topt = 17, Tmax =27)
#----------------------------------------------
# Growth simulation
#---------------------------------------------
# von Bertalaffy : L = f(age)
vonBertalanffyGrowth = function(age, L0, Linf, K){
t0 = log(1 - L0 / Linf) / K
return(Linf * (1 - exp(-K * (age - t0))))
}
# von Bertalanffy inverse : age = f(L)
# with L0 rather the usual t0 parameter
vonBertalanffyInverse = function(L, L0, Linf, K){
t0 = log(1 - L0/Linf)/K
return(t0 - log(1 - L/Linf)/K)
}
# von Bertalanffy increment
# pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche
vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26)
L = matrix(nrow = nStep + 1)
L[1] = L0
for (i in 1:nStep) {
mu = log((Linf - L[i]) * (1 - exp(-K * timeStepDuration))) - sigma * sigma / 2
increment = exp(rnorm(1, mu, sigma))
if (withTempEffect) {
increment = increment * tempEffect[((i - 1) %% 4) + 1]
}
L[i + 1] = L[i] + increment
}
return(L)
}
vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){
if (sigma == 0) {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
increment = exp(mu)
}
else {
mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
increment = exp(rnorm(length(mu), mu, sigma))
}
L = L + increment * tempEffect
return(L)
}
# -------------------------------------------------------
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
# Dispersal
# see see (Chapman et al., 2007; Holloway et al., 2016)
# -------------------------------------------------------
logitKernel = function(distance, alpha0, alpha1, meanInterDistance, standardDeviationInterDistance,
basinSurface = 0, alpha2 = 0, meanBasinSurface = 0, standardDeviationBasinSurface = 1,
fishLength = 0, alpha3 = 0, meanLengthAtRepro = 0, standardDeviationLengthAtRepro = 1) {
logitW = alpha0 - alpha1 * (distance - meanInterDistance) / standardDeviationInterDistance +
alpha2 * (basinSurface - meanBasinSurface) / standardDeviationBasinSurface +
alpha3 * (fishLength - meanLengthAtRepro) / standardDeviationInterDistance;
W = 1 / (1+ exp(-logitW))
return(W)
}
# extended negative exponential kernel
eneKernel = function(distance, alpha_D, beta_D){
W = exp(-alpha_D * distance ^ beta_D)
return(W)
}
strayerSurvival = function(distance, m ){
# m: extra mortality rate for strayers (L-1)
return(exp(-m * distance))
}
#-------------------------------------------------------
# Spwaner survival before reproduction
# Dome-shape curve with temperature effect
#-------------------------------------------------------
spawnerSurvivalPreReproduction <- function(Triver,Tmin, Topt, Tmax, parSurv){
#River survival
survProbOptRiver = parSurv$survProbOptRiver
spRiver = survProbOptRiver * temperatureEffect(Triver, Tmin, Topt, Tmax)
#StockAfterSurv = S * SpRiver
return(spRiver)
}
#----------------------------------------------
# Stock Recruitement relationship
#See Rougier et al, 2014; 2015 and Jatteau et al., 2017
#----------------------------------------------
#survival of larvae 14dph from Jatteau et al, 2017
stockRecruitementRelationship <- function(Triver,Tmin, Topt, Tmax, survivalStock, parRep, parSurv, surfaceWatershed = 84810) {
#surfaceWatershed = 84810 Surface de la Garonne
lambda = parRep$lambda
deltaT = parRep$deltaT
survOptRep = parRep$survOptRep
eta = parRep$eta #simule l'effet allee
ratioS95_S50 = parRep$ratioS95_S50
alpha = parRep$a #species fecundity
survProbOptRiver = parSurv$survProbOptRiver
periodAtSea = 5 - deltaT
###--------------------- SR relationship -----------------
#parametre c de la RS de BH int?grant un effet du BV consid?r?
cj = lambda/surfaceWatershed
#parametre b repr?sentant la mortalit? densit? d?pendante de la RS de BH int?grant un effet de la temperature
# bj = (-(1/deltaTrecruitement))*
# log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0))
bj = -log(survOptRep * temperatureEffect(Triver, Tmin , Topt , Tmax)) / deltaT
#parametre a (f?condit? de l'esp?ce) de la RS de BH int?grant un effet de la temperature
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
alphaj = (bj * exp(-bj * deltaT)) / (cj * (1 - exp(-bj * deltaT)))
alphaj[is.na(alphaj)] <- 0
#Bj = param?tre de la relation SR int?grant l'effet de la temp?rature
betaj = bj/(alpha*cj*(1 - exp(-bj*deltaT)))
#p = proportion de g?niteurs participant ? la reproduction en focntion de la quantit? de g?niteur total
#p = 1/(1+exp(-log(19)*(S-n)/(Ratio*surfaceWatershed)))
S95 = eta * surfaceWatershed
S50 = S95/ratioS95_S50
p = 1/(1 + exp(-log(19)*(survivalStock - S50)/(S95 - S50)))
#relation Stock Recrutement ie calcul le nombre de recrues en fonction du nombre de g?niteurs et de la T en int?grant l'effet Allee
#R0 = aj * S * p
alleeEffect = 1/(1 + exp(-log(19)*(survivalStock - eta/ratioS95_S50*surfaceWatershed)/(eta*surfaceWatershed - eta/ratioS95_S50*surfaceWatershed)))
Rj = (alphaj * survivalStock * alleeEffect)/(betaj + survivalStock * alleeEffect)
#Rj = ((aj * S) * p)/(Bj +S * p)
stockRecruitement = as.vector(Rj)
return(stockRecruitement)
}
# ----------------------------------------------
# Spawner survival after reproduction
#----------------------------------------------
#Logit with temperature effect depending on a Triver
spawnerSurvivalPostReproductionTempRiver <- function(Triver, coeffa, coeffb){
spRiverPostSpawn = (coeffb/ (1 + exp(coeffa*(Triver))))
return(spRiverPostSpawn)
}
#Logit with temperature effect depending on the difference between the river temperature and Tref
spawnerSurvivalPostReproductionTempRef <- function(Triver,Tref, coeffa, coeffb){
spRiverPostSpawn = (coeffb/ (1 + exp(coeffa*(Triver-Tref))))
return(spRiverPostSpawn)
}
#Dome-shape curve with temperature effect
spawnerSurvivalPostReproductionWithBellCurve <- function(Triver, Tmin, Topt, Tmax, coeffa, coeffb){
#P1:
#SpRiverPostSpawn = probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit*(TemperatureEffect(Triver, Tmin, Topt, Tmax))))
#P2:
#SpRiverPostSpawn = probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit))
#SpRiverPostSpawnTempEffect = (probSurvAfterRepro/(probSurvAfterRepro + (1 - probSurvAfterRepro)*exp(-coeffLogit))) * TemperatureEffect(Triver, Tmin, Topt, Tmax)
#P3:
spRiverPostSpawn = coeffb * temperatureEffect(Triver, Tmin, Topt, Tmax)
#stockAfterSpawn = S * spRiverPostSpawn
return(spRiverPostSpawn)
}
211212213