Life history parameters

Get fishbase data

load(url("https://github.com//fishnets//fishnets//blob//master//data//fishbase-web//fishbase-web.RData?raw=True"))

Select turbot as a case study

lh=subset(fb,species=="Psetta maxima")

names(lh)[c(14,17)] = c("l50","a50")
lh=lh[,c("linf","k","t0","a","b","a50","l50")]

head(lh)

     linf    k    t0      a    b a50  l50
3697 53.9 0.27 -0.12     NA   NA  NA   NA
3698 81.5 0.21 -0.48     NA   NA  NA   NA
3699 54.6 0.31 -0.12 0.0189 2.93   4 49.0
3700 78.6 0.19 -0.80     NA   NA  NA   NA
3701 54.4 0.23 -0.22 0.0219 2.92  NA 37.5
3702 73.6 0.28  0.08 0.0189 2.93   4 49.0

Get the means and create an FLPar object

lh=apply(lh,2,mean,na.rm=T)
lh=FLPar(lh)

Values can be replace with any better estimates that may be available.

lhPar fills in missing values using life history theory, while quantities like selection-at-age are set using defaults, in this case set to be the same as maturity-at-age.

par=lhPar(lh)
par

An object of class "FLPar"
params
     linf         k        t0         a         b     ato95       a50 
  58.6840    0.2868   -0.3120    0.0204    2.9250    1.0000    4.0000 
     asym        bg        a1        sl        sr         s         v 
   1.0000    2.9250    5.0000    1.0000 5000.0000    0.9000 1000.0000 
      l50 
  43.2500 
units:  cm 

The parameters are then used by lhEql to simulate the equilibrium dynamics by combining the spawner/yield per recruit relationships with a stock recruitment relationship, by creating an FLBRP object

eq=lhEql(par)

Figure 1 Vectors of m, selection pattern, maturity and weight-at-age.

plot(eq)

Figure 2 Expected, equilibrium, dynamics and reference points.

To model time series the FLBRP object created by lhEql is coerced to an FLStock object and then projected forward

Simulate a time series that represents a stock that was originally lightly exploited, then effort is increased until the stock is overfished and then fishing pressure was reduced to recover the stock biomass to BMSY.

fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(.1,19),
                                              seq(.1,2,length.out = 30)[-30],
                                              seq(2,1.0,length.out = 10),
                                              rep(1.0,61)))[,1:105]
plot(fbar(eq))

Coercing the FLBRP object into an FLStock

om=as(eq,"FLStock")

Call fwd() with the stock, the target F and a stock recruitment relationship

om=fwd(om,fbar=fbar(om)[,-1], sr=eq)

Figure 3 Stochastic Time series of F, SSB, recruitment and yield

plot(FLQuants(om,   
          "ssb" = function(x) ssb(x)%/%refpts( eq)["msy","ssb"], 
          "f" =   function(x) fbar(x)%/%refpts(eq)["msy","harvest"], 
          "rec" = function(x) rec(x)%/%refpts( eq)["msy","rec"], 
          "catch"=function(x) landings(x)%/%refpts(eq)["msy","yield"])) + 
  geom_hline(aes(yintercept=1),col="red") 

Figure 4 Time series relative to MSY benchmarks.

Stochasticity

Add recruitment variability into add the residuals argument, e.g. for 100 Monte Carlo simulations with a CV of 0.3

nits=100

srDev=rlnoise(nits, rec(om) %=% 0, 0.3)

om=propagate(om,nits)

om=fwd(om,fbar=fbar(om)[,-1], sr=eq, residuals=srDev)

plot(om)

Figure 5 Simulate a stock with increasing F