Goals of this practical

By the end of this practical, you should be able to:

Estimate natural mortality (M) using several empirical equations
Compare how different equations can give different answers (uncertainty)
Compute an age-dependent M(a) using a weight-based method (Lorenzen)
Estimate total mortality (Z) using a catch curve
Estimate fishing mortality (F) using Z = M + F

Important: empirical formulas are not “truth”. They are tools.
Different methods often disagree — and that disagreement matters.

0) Load packages

library(ggplot2)
library(dplyr)
library(tidyr)
library(tibble)

1) Life-history inputs (toy example)

Choose values that make biological sense for your species and units.

par <- list(
  Linf = 380,   # VB asymptotic length (cm)
  K    = 0.15,  # VB growth coefficient (1/year)
  t0   = -1,    # VB intercept parameter (years)
  Tmax = 30,    # maximum observed age (years)
  T    = 18,    # mean temperature (°C) used by Pauly (1980)

  # Length-weight relationship: W = c * L^d
  # If L is in cm, W will be in kg for this parameterization.
  cLW = 5.4e-06,
  dLW = 3.14
)

2) Natural mortality (M) — age-invariant (constant M)

2.1 Pauly (1980)

\[ \log_{10}(M) = -0.0066 - 0.279\log_{10}(L_{\infty}) + 0.6543\log_{10}(K) + 0.4634\log_{10}(T) \]

Linf in cm
K in year⁻¹
T in °C

M_pauly1980 <- function(Linf, K, T) {
  10^(-0.0066 - 0.279 * log10(Linf) + 0.6543 * log10(K) + 0.4634 * log10(T))
}

2.2 Hoenig (1983)

A classic longevity relationship (originally for Z; often used as a rough proxy for M):

\[ \ln(Z) = 1.46 - 1.01\ln(T_{max}) \]

M_hoenig1983_fish <- function(Tmax) {
  exp(1.46 - 1.01 * log(Tmax))
}

2.3 Jensen “rules of thumb”

Two common approximate relationships:

\(M = 1.5K\)
\(M = 1.65 / t_{mat}\)

M_jensen_K <- function(K) {
  1.5 * K
}

M_jensen_tmat <- function(tmat) {
  1.65 / tmat
}

Quick maturity guess

If you don’t have maturity age \(t_{mat}\), a very rough assumption is:

\[ t_{mat} \approx T_{max}/3 \]

tmat_guess <- par$Tmax / 3
tmat_guess

## [1] 10

2.4 Calculate M values and survival

Annual survival is:

\[ S = e^{-M} \]

M_table <- tibble(
  Method = c(
    "Pauly 1980 (Linf, K, T)",
    "Hoenig 1983 (longevity)",
    "Jensen (1.5*K)",
    "Jensen (1.65/tmat) (tmat guess)"
  ),
  M = c(
    M_pauly1980(par$Linf, par$K, par$T),
    M_hoenig1983_fish(par$Tmax),
    M_jensen_K(par$K),
    M_jensen_tmat(tmat_guess)
  )
) %>%
  mutate(M_annual_survival = exp(-M))

M_table

## # A tibble: 4 × 3
##   Method                              M M_annual_survival
##   <chr>                           <dbl>             <dbl>
## 1 Pauly 1980 (Linf, K, T)         0.207             0.813
## 2 Hoenig 1983 (longevity)         0.139             0.870
## 3 Jensen (1.5*K)                  0.225             0.799
## 4 Jensen (1.65/tmat) (tmat guess) 0.165             0.848

Plot: compare constant M estimates

ggplot(M_table, aes(x = reorder(Method, M), y = M, fill = Method)) +
  geom_col() +
  coord_flip() +
  theme_bw() +
  labs(
    title = "Natural mortality (M) from different empirical equations",
    x = "Estimator",
    y = "M (year⁻¹)"
  ) +
  theme(legend.position = "none")

3) Natural mortality (M) — age-dependent example (Lorenzen)

3.1 Idea

Small individuals typically have higher natural mortality than large individuals.

One common weight-based relationship is:

\[ M(W) = 3.0 \cdot W^{-0.288} \]

⚠️ Units matter a LOT.
This is often applied with W in grams.

3.2 Length-at-age (VB) and Weight-at-age (L–W)

vb_len <- function(age, Linf, K, t0) {
  Linf * (1 - exp(-K * (age - t0)))
}

lw_weight <- function(L, cLW, dLW) {
  cLW * (L^dLW)
}

age_vec  <- 0:par$Tmax
L_at_age <- vb_len(age_vec, par$Linf, par$K, par$t0)
W_at_age <- lw_weight(L_at_age, par$cLW, par$dLW)  # kg (for this parameterization)

3.3 Lorenzen M(W)

Convert kg → g for the Lorenzen formula:

M_lorenzen <- function(W_g) {
  3.0 * (W_g^(-0.288))
}

W_at_age_g <- W_at_age * 1000  # kg -> g
M_age      <- M_lorenzen(W_at_age_g)

mort_df <- tibble(
  Age    = age_vec,
  Length = L_at_age,
  Weight = W_at_age,
  M_age  = M_age
)

mort_df %>% head(10)

## # A tibble: 10 × 4
##      Age Length Weight  M_age
##    <int>  <dbl>  <dbl>  <dbl>
##  1     0   52.9   1.40 0.373 
##  2     1   98.5   9.81 0.213 
##  3     2  138.   28.1  0.157 
##  4     3  171.   55.9  0.129 
##  5     4  201.   91.4  0.112 
##  6     5  226.  132.   0.101 
##  7     6  247.  176.   0.0926
##  8     7  266.  221.   0.0867
##  9     8  281.  265.   0.0822
## 10     9  295.  308.   0.0788

Plot: M(age) from Lorenzen

ggplot(mort_df, aes(Age, M_age)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2) +
  theme_bw() +
  labs(
    title = "Age-dependent M (Lorenzen weight-based example)",
    x = "Age (years)",
    y = "M(age) (year⁻¹)"
  )

4) Compare M across ages (constant vs age-dependent)

Constant M methods are horizontal lines (same M for every age).
Lorenzen varies with age.

# Expand constant M methods across ages
M_const_by_age <- M_table %>%
  select(Method, M) %>%
  expand_grid(Age = mort_df$Age)

# Lorenzen curve
M_lorenzen_by_age <- mort_df %>%
  transmute(Age = Age,
            Method = "Lorenzen (age-dependent)",
            M = M_age)

# Combine
M_all_age <- bind_rows(M_const_by_age, M_lorenzen_by_age)

ggplot(M_all_age, aes(x = Age, y = M, colour = Method)) +
  geom_line(linewidth = 1.1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  guides(colour = guide_legend(nrow = 2, byrow = TRUE))+
  labs(
    title = "Comparison of natural mortality (M) across ages",
    x = "Age (years)",
    y = "M (year⁻¹)",
    colour = "Method"
  )

5) Compare annual survival across ages

Annual survival:

\[ S = e^{-M} \]

S_all_age <- M_all_age %>%
  mutate(S = exp(-M))

ggplot(S_all_age, aes(x = Age, y = S, colour = Method)) +
  geom_line(linewidth = 1.1) +
  theme_bw() +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    title = "Comparison of annual survival across ages",
    x = "Age (years)",
    y = "Annual survival rate (S)",
    colour = "Method"
  )

6) Catch curve: estimate total mortality Z

6.1 Concept

If recruitment and selectivity are stable and ages are fully recruited:

\[ N(a)\approx e^{-Za} \Rightarrow \log(N(a)) = \text{constant} - Za \]

So the slope of log(N) vs Age is −Z.

6.2 Simulated numbers-at-age dataset

catch_curve <- tibble(
  Age = 1:14,
  N   = c(250, 400, 1200, 1600, 1400, 1100, 900, 650, 450, 310, 220, 160, 120, 80)
)
catch_curve

## # A tibble: 14 × 2
##      Age     N
##    <int> <dbl>
##  1     1   250
##  2     2   400
##  3     3  1200
##  4     4  1600
##  5     5  1400
##  6     6  1100
##  7     7   900
##  8     8   650
##  9     9   450
## 10    10   310
## 11    11   220
## 12    12   160
## 13    13   120
## 14    14    80

Plot: catch-at-age

ggplot(catch_curve, aes(Age, N)) +
  geom_col() +
  theme_bw() +
  labs(
    title = "Catch-at-age (toy example)",
    x = "Age",
    y = "Count (N)"
  )

6.3 Fit a catch curve (choose fully recruited ages)

Choose a cutoff age (students discuss): Age ≥ cut_age

cut_age <- 4
cc_fit <- catch_curve %>% filter(Age >= cut_age)

fit_cc <- lm(log(N) ~ Age, data = cc_fit)
summary(fit_cc)

## 
## Call:
## lm(formula = log(N) ~ Age, data = cc_fit)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.209711 -0.030837 -0.006115  0.055189  0.146213 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.82919    0.09688   91.13 1.17e-14 ***
## Age         -0.31043    0.01016  -30.57 2.10e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1065 on 9 degrees of freedom
## Multiple R-squared:  0.9905, Adjusted R-squared:  0.9894 
## F-statistic: 934.2 on 1 and 9 DF,  p-value: 2.103e-10

Z_hat <- -coef(fit_cc)[["Age"]]
Z_hat

## [1] 0.3104292

Plot: ALL ages + fitted descending line (only fitted ages)

cc_fit <- cc_fit %>%
  mutate(logN_hat = predict(fit_cc, newdata = cc_fit))

ggplot() +
  geom_point(data = catch_curve, aes(Age, log(N)),
             size = 2, alpha = 0.7) +
  geom_line(data = cc_fit, aes(Age, logN_hat),
            linewidth = 1.2) +
  geom_point(data = cc_fit, aes(Age, log(N)),
             colour = "red", size = 2.5) +
  theme_bw() +
  labs(
    title = paste0("Catch curve estimate: Z = ", round(Z_hat, 3), " year⁻¹"),
    subtitle = paste0("All ages shown; model fit only for ages ≥ ", cut_age),
    x = "Age",
    y = "log(N)"
  )

7) Fishing mortality: F = Z − M

7.1 Example using Pauly M

M_pauly <- M_table$M[M_table$Method == "Pauly 1980 (Linf, K, T)"]
M_pauly

## [1] 0.2071334

F_hat <- Z_hat - M_pauly
F_hat

## [1] 0.1032958

7.2 Exercises:

1) Compute and F using ALL constant-M methods (F = Z - M, for each method)

2) Make a plot to compare F across methods

If F becomes negative, it usually means: M is too high relative to Z, or
catch curve assumptions are violated (not fully recruited ages, selectivity changes, variable recruitment, etc.).

8) Age-dependent fishing mortality: F(age) = Z − M(age)

If M is age-dependent, then F can also vary with age:

mort_df <- mort_df %>%
  mutate(
    Z = Z_hat,
    F_age = Z_hat - M_age
  )

mort_df %>% select(Age, M_age, Z, F_age) %>% head(10)

## # A tibble: 10 × 4
##      Age  M_age     Z   F_age
##    <int>  <dbl> <dbl>   <dbl>
##  1     0 0.373  0.310 -0.0623
##  2     1 0.213  0.310  0.0978
##  3     2 0.157  0.310  0.153 
##  4     3 0.129  0.310  0.182 
##  5     4 0.112  0.310  0.199 
##  6     5 0.101  0.310  0.210 
##  7     6 0.0926 0.310  0.218 
##  8     7 0.0867 0.310  0.224 
##  9     8 0.0822 0.310  0.228 
## 10     9 0.0788 0.310  0.232

Plot: M(age) and F(age)

ggplot(mort_df, aes(x = Age)) +
  geom_line(aes(y = M_age, colour = "M(age)"), linewidth = 1.2) +
  geom_line(aes(y = F_age, colour = "F(age) = Z - M(age)"), linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  labs(
    title = "Age-dependent natural mortality and fishing mortality",
    subtitle = paste0("Using catch-curve Z = ", round(Z_hat, 3)),
    x = "Age (years)",
    y = "Instantaneous rate (year⁻¹)",
    colour = "Curve"
  )

Interpretation note: F(age) is only meaningful for ages that are fully recruited to the fishery
(in this toy example: ages ≥ cut_age).

Wrap-up

Key message: M is uncertain and very influential in stock assessments. Uncertainty propagates to F.
Good assessments explore sensitivity to M assumptions.

cat("\nCONGRATS!\n")

## 
## CONGRATS!

TP3 (Part 2): Mortality (M, Z and F) from Life-History

Rui Coelho (rpcoelho@ualg.pt)

Jan 2026