Fish biomass production modeling
Matheus de Barros
2023-04-05
Introduction
Accumulated fish biomass production over time can be determined by applying known life history relationships such as length-weight, mortality, and growth. This framework is based on modified equations by zu Ermgassen et al. (2015).
First, the rate of production can be determined by predicting the decline in numbers of a cohort over time and the predicted biomass for an individual fish at a given time:
\[ \displaystyle \frac{\partial P}{\partial t} = N(t)\frac{\partial W}{\partial t} \] \[ \frac{\partial N}{\partial t} = e^{-M(t)} \] where P symbolizes production for a whole population or cohort, N(t) is numbers at a given time given the rate of natural mortality (M), and W is predicted individual weight at a given time.
Individual weight can be predicted by applying coefficients of a length-weight relationship to the predicted size over time:
\[ \frac{\partial W}{\partial t} = \begin{cases} \ a(L_{t-1}+GR)^{b},& \text{if } L < L_{m}\\ a[L_{inf} e^{-k(t-t_{0})}]^{b}, \text{otherwise} \end{cases} \ \]
in this equation, growth follows a linear trajectory over time (represented by mm/day, GR) until the fish reaches its age at maturation (Lm). After maturation, growth is predicted by the classical von Bertalanffy equation. Coefficients of the length-weight relationship (a and b) are estimated previously.
When the first equation is integrated over time, it allows us to estimate gross biomass production (GP) over a given period. This accounts for every biomass produced in the specified time frame, including live and dead individuals:
\[ GP = \int_{t_{r}}^{t_{max}}\frac{\partial P}{\partial t} \]
\[ GP = \int_{t_{r}}^{t_{max}}\frac{\partial W}{\partial t}\frac{\partial N}{\partial t} \] I’ve written the following R code to estimate gross production over time using the above framework:
#-------------------------------------------------------------------------------
# FUNCTIONS TO MODEL BIOMASS ACCUMULATION
#-------------------------------------------------------------------------------
set.seed(2023)
pred_B <- function(GR = NULL, N0 = NULL, M_ref = NULL, L_ref = NULL,
Lm = NULL, K = NULL, Linf = NULL, t0 = NULL,
a = NULL, b = NULL, L0 = NULL, K_sup = NULL,
t.steps = NULL, plot = FALSE) {
# growth is assumed to follow a linear trajectory at the juvenile phase
# and a von Bertalanffy type after maturation
L <- rep(NA, t.steps)
L_prior <- rep(NA, t.steps)
for(t in 2:length(L)) {
L_prior[1] <- L0
L_prior[t] <- L_prior[t-1]+GR
}
for (t in 1:length(L)) {
L[t] <- ifelse(L_prior[t] < Lm,
L_prior[t]*1,
(Linf*(1-exp(- (K/365) *(t- (t0/365))))))
}
point <- c(which.min(abs(L - Lm))) # find value closest to breakpoint Lm
L_corr <- rep(NA, t.steps)
for(t in 1:length(L_corr)) {
#apply a correction factor to make VB growth continue from the last point
L_corr[t] <- ifelse(L_prior[t] < Lm,
L_prior[t],
(Linf*(1-exp(-(K/365)*(t-(t0/365))))+abs(L_corr[point]-L[point+1])))
}
L_VB <- rep(NA, t.steps)
for(t in 1:length(L_VB)) {
L_VB[t] <- Linf*(1-exp(-K*(t)))
}
# estimation of mortality
# mortality is size-dependent and is calculated based on known values for a given length
M <- rep(NA, t.steps)
for(t in 1:length(M)) {
M[t] <- (M_ref/365)*(L_ref/L_corr[t])
}
# estimation of numbers at age
N <- rep(NA, t.steps)
for(t in 2:length(N)) {
N[1] <- N0
N[t] <- N[t-1]*exp(-M[t])
}
# calculate biomass at size
B <- rep(NA, t.steps)
for(t in 1:length(B)) {
B[t] <- (a*L_corr[t])^b
}
# calculate cohort biomass
CB <-rep(NA, t.steps)
for(t in 1:length(CB)) {
CB[t] <- N[t] * B[t]
}
P_acc <- cumsum(CB) # derived quantities
P_accK <- rep(NA, t.steps)
for(t in 1:length(P_accK)) {
P_accK[t] <- P_acc[t]*(1-P_acc[t]/K_sup)# some code to include a carrying capacity.
# Yet to be developed
}
if (plot == TRUE) {
plot(P_acc, type = "l")
}
output <- data.frame(L_corr, N, B, CB, P_acc)
return(output)
}Observe that natural mortality in this function is modeled with a size-dependent relationship from Dureuil et al. (2021), so every size at a given time (Lt) is going to have its own mortality. It needs a reference value (Mr) at a given size (Lr):
\[ M(t) = M_{r}(L_{r}/L_{t}) \]
Let’s take a look at the output with a simple example. Let growth rate (GR) be = 0.5 mm/day, reference mortality rate for a 120 mm fish be 1.6 day-1, length at maturity = 141 mm, von Bertalanffy growth rate (k) = 0.33 day-1, asymptotic length (Linf) = 200 mm, and initial size (L0) = 30 mm. Initial abundance at t = 0 (N0) is 0.5 individuals/m2, so output in biomass is going to be in g/m2/year.
out <- pred_B(GR = 0.5, M_ref = 1.6, L_ref = 120, Lm = 141, K = 0.33/365,
Linf = 200, t0 = -1.1, t.steps = 365*5, L0 = 30, N0 = 0.5,
a = 0.0064, b = 3.2, plot = TRUE, K_sup = 10)
In this case, gross fish production asymptotes at around 25g of fish/m2. Let’s decrease the growth rate at the juvenile phase by 0.1 mm/day and see what happens:
out2 <- pred_B(GR = 0.4, M_ref = 1.6, L_ref = 120, Lm = 141, K = 0.33/365,
Linf = 200, t0 = -1.1, t.steps = 365*5, L0 = 30, N0 = 0.5,
a = 0.0064, b = 3.2, plot = TRUE, K_sup = 10)
The next step should be to obtain measures of variability on the output. For that, I’ve written a Monte Carlo framework to incorporate variation in input parameters in the production calculations:
#-------------------------------------------------------------------------------
# RUN A MONTE CARLO SIMULATION TO ESTIMATE UNCERTAINTY
#-------------------------------------------------------------------------------
stoch_P <- function(GR.mu = NULL, GR.sd = NULL,
N0.mu = NULL, N0.sd = NULL,
M_ref.mu = NULL, M_ref.sd = NULL,
L_ref.mu = NULL, L_ref.sd = NULL,
Lm.mu = NULL, Lm.sd = NULL,
K.mu = NULL, K.sd = NULL,
Linf.mu = NULL, Linf.sd = NULL,
t0.mu = NULL, t0.sd = NULL,
lwa.mu = NULL, lwa.sd = NULL,
lwb.mu = NULL, lwb.sd = NULL,
L0.mu = NULL, L0.sd = NULL,
K_sup.mu = NULL, K_sup.sd = NULL,
t.steps = NULL, N = 100,
plot = FALSE, N_years = NULL) {
# create a data frame with distributions of input parameters to be looped over
prior.df <- data.frame(rnorm(N, GR.mu, GR.sd),
rnorm(N, N0.mu, N0.sd),
rnorm(N, M_ref.mu, M_ref.sd),
rnorm(N, L_ref.mu, L_ref.sd),
rnorm(N, Lm.mu, Lm.sd),
rnorm(N, K.mu, K.sd),
rnorm(N, Linf.mu, Linf.sd),
rnorm(N, t0.mu, t0.sd),
rnorm(N, lwa.mu, lwa.sd),
rnorm(N, lwb.mu, lwb.sd),
rnorm(N, L0.mu, L0.sd),
rnorm(N, K_sup.mu, K_sup.sd))
names(prior.df) <- c('GR', 'N0', 'M_ref', 'L_ref', 'Lm', 'K', 'Linf', 't0', 'lwa', 'lwb', 'L0', 'K_sup')
out.df <- NULL
pb <- txtProgressBar(min = 0, max = N, style = 3, char = "-")
for(i in 1:N) {
out <- pred_B(GR = prior.df$GR[i],
N0 = prior.df$N0[i],
M_ref = prior.df$M_ref[i],
L_ref = prior.df$L_ref[i],
Lm = prior.df$Lm[i],
K = prior.df$K[i],
Linf = prior.df$Linf[i],
t0 = prior.df$t0[i],
a = prior.df$lwa[i],
b = prior.df$lwb[i],
L0 = prior.df$L0[i],
K_sup = prior.df$K_sup[i], t.steps = 365*N_years)
out.df <- cbind(out.df, out$P_acc, deparse.level = 0)
setTxtProgressBar(pb, i)
}
Bacc.mu <- apply(out.df, 1, median, na.rm = TRUE)
Bacc.up <- apply(out.df,1,quantile, probs = c(.9), na.rm = TRUE)
Bacc.low <- apply(out.df,1,quantile, probs = c(.1), na.rm = TRUE)
plot(Bacc.mu, type = "l", lty = 1, lwd = 1.5, xlab = "Time(days)", ylab = "GP(g/m2)")
lines(Bacc.low, lty = 3)
lines(Bacc.up, lty = 3)
Bacc.out <- data.frame(Bacc.mu, Bacc.low, Bacc.up)
close(pb)
return(Bacc.out)
}Now let’s test it. The number of iterations in this example is set to N = 100, but a much higher value should be set for when using real data so the function can loop over the whole distributions of input parameters.
out3 <- stoch_P(GR.mu = 0.4, GR.sd = 0.001, N0.mu = 0.3, N0.sd = 0.001,
M_ref.mu = 1.6, M_ref.sd = 0.1, L_ref.mu = 120, L_ref.sd = 0.1,
Lm.mu = 120, Lm.sd = 10, K.mu = 0.33, K.sd = 0.005,
Linf.mu = 200, Linf.sd = 10, t0.mu = 0, t0.sd = 0.01,
lwa.mu = 0.0064, lwa.sd = 0.0000001, lwb.mu = 3.2,
lwb.sd = 0.001, L0.mu = 30, L0.sd = 0.1, K_sup.mu = 40,
K_sup.sd = 0, N = 100, N_years = 7)##
|
| | 0%
|
|- | 1%
|
|- | 2%
|
|-- | 3%
|
|--- | 4%
|
|---- | 5%
|
|---- | 6%
|
|----- | 7%
|
|------ | 8%
|
|------ | 9%
|
|------- | 10%
|
|-------- | 11%
|
|-------- | 12%
|
|--------- | 13%
|
|---------- | 14%
|
|---------- | 15%
|
|----------- | 16%
|
|------------ | 17%
|
|------------- | 18%
|
|------------- | 19%
|
|-------------- | 20%
|
|--------------- | 21%
|
|--------------- | 22%
|
|---------------- | 23%
|
|----------------- | 24%
|
|------------------ | 25%
|
|------------------ | 26%
|
|------------------- | 27%
|
|-------------------- | 28%
|
|-------------------- | 29%
|
|--------------------- | 30%
|
|---------------------- | 31%
|
|---------------------- | 32%
|
|----------------------- | 33%
|
|------------------------ | 34%
|
|------------------------ | 35%
|
|------------------------- | 36%
|
|-------------------------- | 37%
|
|--------------------------- | 38%
|
|--------------------------- | 39%
|
|---------------------------- | 40%
|
|----------------------------- | 41%
|
|----------------------------- | 42%
|
|------------------------------ | 43%
|
|------------------------------- | 44%
|
|-------------------------------- | 45%
|
|-------------------------------- | 46%
|
|--------------------------------- | 47%
|
|---------------------------------- | 48%
|
|---------------------------------- | 49%
|
|----------------------------------- | 50%
|
|------------------------------------ | 51%
|
|------------------------------------ | 52%
|
|------------------------------------- | 53%
|
|-------------------------------------- | 54%
|
|-------------------------------------- | 55%
|
|--------------------------------------- | 56%
|
|---------------------------------------- | 57%
|
|----------------------------------------- | 58%
|
|----------------------------------------- | 59%
|
|------------------------------------------ | 60%
|
|------------------------------------------- | 61%
|
|------------------------------------------- | 62%
|
|-------------------------------------------- | 63%
|
|--------------------------------------------- | 64%
|
|---------------------------------------------- | 65%
|
|---------------------------------------------- | 66%
|
|----------------------------------------------- | 67%
|
|------------------------------------------------ | 68%
|
|------------------------------------------------ | 69%
|
|------------------------------------------------- | 70%
|
|-------------------------------------------------- | 71%
|
|-------------------------------------------------- | 72%
|
|--------------------------------------------------- | 73%
|
|---------------------------------------------------- | 74%
|
|---------------------------------------------------- | 75%
|
|----------------------------------------------------- | 76%
|
|------------------------------------------------------ | 77%
|
|------------------------------------------------------- | 78%
|
|------------------------------------------------------- | 79%
|
|-------------------------------------------------------- | 80%
|
|--------------------------------------------------------- | 81%
|
|--------------------------------------------------------- | 82%
|
|---------------------------------------------------------- | 83%
|
|----------------------------------------------------------- | 84%
|
|------------------------------------------------------------ | 85%
|
|------------------------------------------------------------ | 86%
|
|------------------------------------------------------------- | 87%
|
|-------------------------------------------------------------- | 88%
|
|-------------------------------------------------------------- | 89%
|
|--------------------------------------------------------------- | 90%
|
|---------------------------------------------------------------- | 91%
|
|---------------------------------------------------------------- | 92%
|
|----------------------------------------------------------------- | 93%
|
|------------------------------------------------------------------ | 94%
|
|------------------------------------------------------------------ | 95%
|
|------------------------------------------------------------------- | 96%
|
|-------------------------------------------------------------------- | 97%
|
|--------------------------------------------------------------------- | 98%
|
|--------------------------------------------------------------------- | 99%
|
|----------------------------------------------------------------------| 100%
Yay! Now we’ve estimated the gross production for this hypothetical species with uncertainty. The final values in the output can now be used as posterior distributions (mean + 90% C.I.).
References
Dureuil, M., Aeberhard, W. H., Burnett, K. A., Hueter, R. E., Tyminski, J. P., & Worm, B. (2021). Unified natural mortality estimation for teleosts and elasmobranchs. Marine Ecology Progress Series, 667, 113-129.
Zu Ermgassen, P. S., Grabowski, J. H., Gair, J. R., & Powers, S. P. (2016). Quantifying fish and mobile invertebrate production from a threatened nursery habitat. Journal of Applied Ecology, 53(2), 596-606.