GR3Dfunction.R 7.41 KiB
#----------------------------------------------
# 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