GR3Dfunction.R 7.08 KB
Newer Older
1
2
3
#----------------------------------------------
# Temperature effect
#----------------------------------------------
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's avatar
patrick.lambert committed
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)
22

23
24
25
#----------------------------------------------
# Growth simulation
#---------------------------------------------
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's avatar
patrick.lambert committed
41
vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
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's avatar
patrick.lambert committed
46
    mu = log((Linf - L[i]) * (1 - exp(-K * timeStepDuration))) - sigma * sigma / 2
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's avatar
patrick.lambert committed
56
vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){
57
  if (sigma == 0) {
patrick.lambert's avatar
patrick.lambert committed
58
    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
59
60
61
    increment = exp(mu)
  }
  else {
patrick.lambert's avatar
patrick.lambert committed
62
    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
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;
  
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))
}
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)
}