GR3Dfunction.R 7.08 KB
 Poulet Camille committed Feb 25, 2021 1 2 3 ``````#---------------------------------------------- # Temperature effect #---------------------------------------------- `````` patrick.lambert committed Dec 05, 2020 4 5 6 7 8 9 10 11 12 ``````# 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) } `````` patrick.lambert committed Feb 17, 2021 13 14 15 16 17 18 19 20 21 ``````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) `````` patrick.lambert committed Dec 05, 2020 22 `````` `````` Poulet Camille committed Feb 25, 2021 23 24 25 ``````#---------------------------------------------- # Growth simulation #--------------------------------------------- `````` patrick.lambert committed Dec 05, 2020 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 ``````# 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 `````` patrick.lambert committed Feb 17, 2021 41 ``````vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){ `````` patrick.lambert committed Dec 05, 2020 42 43 44 45 `````` tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26) L = matrix(nrow = nStep + 1) L[1] = L0 for (i in 1:nStep) { `````` patrick.lambert committed Feb 17, 2021 46 `````` mu = log((Linf - L[i]) * (1 - exp(-K * timeStepDuration))) - sigma * sigma / 2 `````` patrick.lambert committed Dec 05, 2020 47 48 49 50 51 52 53 54 55 `````` increment = exp(rnorm(1, mu, sigma)) if (withTempEffect) { increment = increment * tempEffect[((i - 1) %% 4) + 1] } L[i + 1] = L[i] + increment } return(L) } `````` patrick.lambert committed Feb 17, 2021 56 ``````vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){ `````` patrick.lambert committed Dec 05, 2020 57 `````` if (sigma == 0) { `````` patrick.lambert committed Feb 17, 2021 58 `````` mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) `````` patrick.lambert committed Dec 05, 2020 59 60 61 `````` increment = exp(mu) } else { `````` patrick.lambert committed Feb 17, 2021 62 `````` mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2 `````` patrick.lambert committed Dec 05, 2020 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 `````` 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; `````` patrick.lambert committed Dec 18, 2020 80 81 82 83 84 85 86 87 88 89 90 91 92 `````` 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)) } `````` Poulet Camille committed Feb 25, 2021 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 `````` #------------------------------------------------------- # 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 intgrant un effet du BV considr cj = lambda/surfaceWatershed #parametre b reprsentant la mortalit densit dpendante de la RS de BH intgrant 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 (fcondit de l'espce) de la RS de BH intgrant un effet de la temperature alphaj = (bj * exp(-bj * deltaT)) / (cj * (1 - exp(-bj * deltaT))) alphaj[is.na(alphaj)] <- 0 #Bj = paramtre de la relation SR intgrant l'effet de la temprature betaj = bj/(alpha*cj*(1 - exp(-bj*deltaT))) #p = proportion de gniteurs participant la reproduction en focntion de la quantit de gniteur 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 gniteurs et de la T en intgrant 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) } ``````