#---------------------------------------------- # 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(1, mu, sigma)) } L = L + increment * tempEffect return(L) } # ------------------------------------------------------- # 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){ #River survival survProbOptRiver = survivePar\$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) { surfaceWatershed = 84810 #Surface de la Garonne lambda = reproducePar\$lambda deltaT = reproducePar\$deltaT survOptRep = reproducePar\$survOptRep n = reproducePar\$eta #simule l'effet allee ratioS95_S50 = reproducePar\$ratioS95_S50 alpha = reproducePar\$a #species fecundity survProbOptRiver = survivePar\$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 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 = n * 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 - n/ratioS95_S50*surfaceWatershed)/(n*surfaceWatershed - n/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 spawnerSurvivalPostReproduction <- 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) }