Introduction

The southern hake (Merluccius australis) is a deepwater demersal species of high commercial value in New Zealand’s fisheries, occurring predominantly south of 40°S within the Exclusive Economic Zone (EEZ). As a member of the Merlucciidae family, hake exhibit biological and ecological traits that make them particularly relevant for stock assessment and sustainable fisheries management. Their life history characteristics—such as moderate longevity (up to 25–30 years), delayed maturity (4–8 years), sexual dimorphism, and distinct growth rates across regions—have direct implications for how these populations are managed and monitored.

In New Zealand, M. australis inhabits continental slope waters primarily between 250 and 900 m in depth, though adults can be found as deep as 1400 m. Juveniles tend to occupy shallower inshore waters (<250 m), where they are vulnerable to trawl fisheries (Horn, 2015). Females grow significantly larger than males, with mature individuals reaching up to 125 cm in total length. Hake are taken almost exclusively by large commercial trawlers and are targeted directly or caught as bycatch in fisheries targeting hoki (Macruronus novaezelandiae) and southern blue whiting (Micromesistius australis). This overlap in fishing effort across species places additional pressure on hake stocks and necessitates careful assessment of stock-specific dynamics.

Hake stock structure in New Zealand waters is divided into three spatially and biologically distinct units, managed independently under the Quota Management System (QMS): the West Coast South Island (HAK 7), the Sub-Antarctic (HAK 1), and the Chatham Rise (HAK 4). HAK 7 historically accounts for the largest annual catches, while the Sub-Antarctic stock (HAK 1) is the most abundant in terms of unfished spawning biomass (94,200 t). Conversely, the Chatham Rise stock (HAK 4) has been the most heavily exploited, with an estimated unfished biomass of just 37,000 t and has been undergoing a rebuilding phase in recent years (Horn, 2015; Fisheries NZ, 2024).

Recent stock assessments suggest that none of the hake stocks currently fall below the management target of 40% of the unfished spawning biomass (B₀), though biological indicators such as maturity-at-age, mortality rates, and growth parameters vary regionally and between sexes. For example, von Bertalanffy growth parameters differ significantly across the three stocks, with females in the Chatham Rise estimated to have K = 0.229 yr⁻¹ and L∞ = 106.5 cm, compared to K = 0.280 yr⁻¹ and L∞ = 99.6 cm for females on the West Coast South Island (Arancibia et al., 2015). These regional differences are critical for constructing accurate population models and informing sustainable catch limits.

From an economic perspective, New Zealand hake fetches a high export price, especially as frozen fillet (up to US$7.73/kg), highlighting its value within the global seafood market (Arancibia et al., 2015). However, the fishery also faces challenges from environmental change. Modelling work by Datta et al. (2024) suggests that deep-sea ecosystems such as the Chatham Rise—home to one of the main hake stocks—may be particularly vulnerable to climate-induced shifts in biomass, especially under warming and fishing scenarios. These findings underscore the need for predictive, data-driven approaches to manage the species in the face of both anthropogenic and environmental pressures.

This report aims to assess the current status and biological parameters of New Zealand hake using a combination of empirical data and modeling techniques. Through the application of standard stock assessment procedures—including the estimation of length-weight relationships, von Bertalanffy growth curves, maturity-at-age, catch-curve mortality analysis, and yield-per-recruit modeling—this report will build a comprehensive understanding of the species’ population dynamics. Model outputs will be used to evaluate current fishing pressure, assess biological sustainability (including spawner-per-recruit), and identify optimal harvest strategies under varying selectivity scenarios.

The following section presents the results of these analyses using R, with annotated code, model fits, residual checks, and interpretations included throughout.

Data Input and Parameter Development

Data Cleaning The raw data has been sourced. The next step is to make all the required factors into numeric variables, and to then remove all NAs from the dataset.

getwd()
## [1] "/Users/willroddick/Desktop/R/410"
hake_data<-read.csv("hake.csv",header = TRUE)
hake_data$age<-as.numeric(hake_data$age)
hake_data$lgth<-as.numeric(hake_data$lgth)
hake_data$weight<-as.numeric(hake_data$weight)
hake_data$sex<-as.factor(hake_data$sex)
hake_data$gonad_stage<-as.numeric(hake_data$gonad_stage)
hake<-subset(hake_data,!is.na(hake_data$age)&!is.na(hake_data$lgth)&!is.na(hake_data$weight)&!is.na(hake_data$sex)
            &!is.na(hake_data$gonad_stage))

Length-Weight Firstly, we craete the length weight subset, and make sure that any sexes other than male/female are removed. Thi is then followed by the first plot.

datLW<-subset(hake,!is.na(hake$lgth)&!is.na(hake$weight))
datLW_sexclean<-subset(datLW,sex%in%c("1","2"))
datLW_sexclean$lgth<-as.numeric(datLW_sexclean$lgth)
datLW_sexclean$weight<-as.numeric(datLW_sexclean$weight)
plot(weight~lgth, datLW_sexclean,pch=19,xlab="Total Length (cm)",ylab="Weight (g)", col="purple", main = "Hake Combined Length-Weight")

Outliers: Immediately we notice 3 present outliers. These 3 outliers are then located and removed from further modelling. This is to ensure the accuracy and stability of the data, and the prevention of distortion for the results. The resulting plot displays no outliers.

datLWbetter1<-datLW_sexclean[-which(datLW_sexclean$weight>15000 & datLW_sexclean$weight<17500 & datLW_sexclean$lgth>85 & datLW_sexclean$lgth<95),]
datLWbetter2<-datLWbetter1[-which(datLWbetter1$weight>2000 & datLWbetter1$weight<3500 & datLWbetter1$lgth>80 & datLWbetter1$lgth<85),]
datLWbest<-datLWbetter2[-which(datLWbetter2$weight>0 & datLWbetter2$weight<2000 & datLWbetter2$lgth>80 & datLWbetter2$lgth<90),]
plot(weight~lgth, datLWbest, pch=19,xlab="Total Length (cm)",ylab="Weight (g)", col="purple", main = "Hake Combined Length-Weight (Outliers Removed)")

Non Linear Squares (NLS) is a method used to find the best-fit parameters for a model that’s not linear in its parameters, by minimizing the sum of squared differences between the model’s predictions and the actual data. Here, we utilize NLS to create a best fit line to insert into the graph, showing a more accurate length-weight relationship for hake. The resulting length-weight plot demonstrates a positive linear relationship, with lengths exceeding 120cm, and weight exceeding 15kg. This aligns with plenary biology information. a = is the scaling factor. b = the exponent of how the length-weight relationship is changed. This typically is ~3 (cubed) as many 3 dimensional animals are assessed using their volume. To assess the level of error present, and whether the model is accurate enough for our graph, a plot displaying the residual and fitted values is created. The resulting plot displays funneling of the data points; this indicates more error occurs as the length-weight increases.

fittedLW <- nls(weight ~ a*lgth^b, data=datLWbest, start=c(a=0.01,b=3))
plot(weight~lgth, datLWbest, pch=19,xlab="Total Length (cm)",ylab="Weight (g)", col="purple", main = "Hake Combined Length-Weight")
lines(predict(fittedLW)~datLWbest$lgth, col="orange",lwd=2)
text(x=40, y=12000, paste("a=", round(coef(fittedLW)[1],3),sep=""))
text(x=45, y=10000, paste("b=", round(coef(fittedLW)[2],3),sep=""))

resLW<-residuals(fittedLW)
fitLW<-fitted(fittedLW)
plot(resLW~fitLW, xlab="Fitted Values",ylab="Residual Values",pch = 19,col="seagreen",main="Residuals Plot Length-Weight")

Males and Females: The same length weight plot, NLS fit and residual error plots were created for both male and female data subsets. The males are notably shorter and lighter, only exceeding 100cm and 10kg respectively. Females are larger, exceeding 120cm and 15kg. These sex-based differences align with the plenary biology statement. There was a potential high outlier observed for the male L-W plot, but was deemed to be acceptable and not in a position to distort. Both male and female plots displayed the same linear relationship and a/b parameters, which gives us confidence in the data. However, both sexes displayed similar residual error funneling trends.

malesLW<-subset(datLWbest,sex =="1")
malesLWfitting <- nls(weight ~ a*lgth^b, data=malesLW, start=c(a=0.01,b=3))
plot(weight~lgth, malesLW, pch=19, xlab="Total Length (cm)",ylab="Weight (g)", col="purple", main = "Hake Male Combined Length-Weight")
lines(predict(malesLWfitting)~malesLW$lgth, col="orange",lwd=2)
text(x=40, y=8500, paste("a=", round(coef(malesLWfitting)[1],3),sep=""))
text(x=45, y=7500, paste("b=", round(coef(malesLWfitting)[2],3),sep=""))

femalesLW<-subset(datLWbest,sex =="2")
femalesLWfitting <- nls(weight ~ a*lgth^b, data=femalesLW, start=c(a=0.01,b=3))
plot(weight~lgth, femalesLW, pch=19, xlab="Total Length (cm)",ylab="Weight (g)",col="purple",main = "Hake Female Combined Length-Weight")
lines(predict(femalesLWfitting)~femalesLW$lgth, col="orange",lwd=2)
text(x=60, y=12500, paste("a=", round(coef(femalesLWfitting)[1],3),sep=""))
text(x=68, y=9000, paste("b=", round(coef(femalesLWfitting)[2],3),sep=""))

resmaleLW<-residuals(malesLWfitting)
fitmaleLW<-fitted(malesLWfitting)
plot(resmaleLW~fitmaleLW, xlab="Fitted Values",ylab="Residual Values",pch = 19,col="seagreen",main="Residuals Plot Length-Weight (Male)")

resfemLW<-residuals(femalesLWfitting)
fitfemaleLW<-fitted(femalesLWfitting)
plot(resfemLW~fitfemaleLW, xlab="Fitted Values",ylab="Residual Values",pch = 19,col="seagreen",main="Residuals Plot Length-Weight (Female)")

Log transformation and Back Transformed Model:

To adjust for the residual error, the combined hake dataset underwent log transformation to attain more accurate parameters and a better fit model by stabilizing variance.

The log-transformed plot demonstrates a clear straight line, as the model has been successfully normalized. We can also see the resulting residual plot has far less funneling than previously observed, indicating that the error has reduced. After proving that the log transformation was successful, we back-transformed the the optimized model. This was then plotted alongside the original model, and the original model utilizing only the unique ages from the hake dataset.

As seen in the plot, all three length-weight models are near identical in curve and position, giving us confidence in the original fitted model. The a & b parameters also remain near identical.

datLWbest$logW<-log10(datLWbest$weight)
datLWbest$logL<-log10(datLWbest$lgth)
plot(logW~logL,data=datLWbest,pch=19,col="darkblue",main="Hake Length-Weight Log-Transformed")
text(x=1.6, y=3.75, paste("a=0.002",sep=""))
text(x=1.65, y=3.4, paste("b=3.266",sep=""))
fitlog<-lm(logW~logL,data=datLWbest)
lines(predict(fitlog)~datLWbest$logL,col="red")

res_log_fit<-residuals(fitlog)
fit_log_fit<-fitted(fitlog)
plot(res_log_fit~fit_log_fit,xlab="Log-Fitted Values",ylab="Log-Residual Values",pch = 19,col="seagreen",main="Log-Residuals Plot Length-Weight")

plot(weight~lgth, datLWbest, pch=19,xlab="Total Length (cm)",ylab="Weight (g)", col="purple3", main = "Hake Combined Length-Weight (Back-Transformed)")
logy<-predict(fitlog)
btx<-10^datLWbest$logL
bty<-10^logy
lines(bty~btx,lwd=2,col="red")
lines(predict(fittedLW)~datLWbest$lgth,col="orange",lwd=2)
xfit<-seq(1:max(datLWbest$lgth))
lines(predict(fittedLW,data.frame(lgth=xfit))~xfit,col="green",lwd=2)
text(x=40, y=11000, paste("a=0.002",sep=""))
text(x=45, y=9000, paste("b=", round(coef(fitlog)[2],3),sep=""))
legend("topleft",
       legend = c("Back-Transformed Model", "Original Model","Original (Unique Ages)"),
       col = c("red", "orange","green"),
       lwd = 2)

Growth

A subset of the data is created using the age and length factors. To create the growth model, we utilize the von Bertalanffy Growth Function (VBGF). The VBGF equation has 3 parameters: Linf = avergae maximum length, K = growth rate and t0 = age at length 0. The VBGF model is then fit to the plot using the NLS function. To further refine and optimize the model, we only use the unique ages (1-26) to prevent repeats of the same age. No outliers were deemed to be distortive, and the parameters were acceptable. To assess the level of error present, and whether the model is accurate enough for our graph, a plot displaying the residual and fitted values is created. The resulting plot displays funneling of the data points; this indicates more error occurs as the length-weight increases. Therefore it is important to interpret the data with caution.

datG1<-subset(hake,!is.na(hake$age)&!is.na(hake$lgth))
datG<-subset(datG1,sex%in%c("1","2"))
vbgf<-nls(lgth~Linf*(1-exp(-K*(age-t0))),data=datG,start=c(Linf=113.65,K=0.163,t0=-1.21))
plot(lgth~age, datG, xlab="Age (Years)",ylab="Total Length (cm)", pch=19, col="69", main="Hake Growth Combined")
text(x=8, y=60, paste("linf = ", round(coef(vbgf)[1],3),sep=""))
text(x=13, y=60, paste("K = ", round(coef(vbgf)[2],3),sep=""))
text(x=18, y=60, paste("t0 = ", round(coef(vbgf)[3],3),sep=""))
uniqueages<-seq(range(datG$age)[1],range(datG$age)[2],1)
lines(predict(vbgf,data.frame(age=uniqueages))~uniqueages,col="red",lwd=2)

res_vb<-residuals(vbgf)
fit_vb<-fitted(vbgf)
plot(res_vb~fit_vb,pch=19,col="seagreen",xlab="Fitted Values",ylab="Residual Values",main="Hake Growth Residuals Plot")

Males & Females The same application of the VBGF model using the NLS function was applied to male and female hake data subsets, along with corresponding residual error plots. The K value (the growth rate) for females was lowered by 0.01 to allow for a more aligned plot. Female hake have higher maximum age and length than male hake, which is in line with prior hake biology. Both male/female residual error plots showed some form of funneling, indicating that the resulting parameters should be interpreted with caution.

malesG<-subset(datG,sex =="1")
vbgfmales<-nls(lgth~Linf*(1-exp(-K*(age-t0))),data=malesG,start=c(Linf=61,K=0.4,t0=0))
plot(lgth~age, malesG, xlab="Age (Years)",ylab="Total Length (cm)", pch=19, col="69", main="Hake Growth Males")
text(x=7, y=60, paste("linf = ", round(coef(vbgfmales)[1],3),sep=""))
text(x=11, y=60, paste("K = ", round(coef(vbgfmales)[2],3),sep=""))
text(x=16, y=60, paste("t0 = ", round(coef(vbgfmales)[3],3),sep=""))
uniqueages_males<-seq(range(malesG$age)[1],range(malesG$age)[2],1)
lines(predict(vbgfmales,data.frame(age=uniqueages_males))~uniqueages_males,col="red",lwd=2)

femalesG<-subset(datG,sex =="2")
plot(lgth~age, femalesG, xlab="Age (Years)",ylab="Total Length (cm)", pch=19, col="69", main="Hake Growth Females")
vbgffemales<-nls(lgth~Linf*(1-exp(-K*(age-t0))),data=femalesG,start=c(Linf=61,K=0.39,t0=0))
text(x=9, y=70, paste("linf = ", round(coef(vbgffemales)[1],3),sep=""))
text(x=14, y=70, paste("K = ", round(coef(vbgffemales)[2],3),sep=""))
text(x=19, y=70, paste("t0 = ", round(coef(vbgffemales)[3],3),sep=""))
#note that we had to change "K" value
uniqueages_females<-seq(range(femalesG$age)[1],range(femalesG$age)[2],1)
lines(predict(vbgffemales,data.frame(age=uniqueages_females))~uniqueages_females,col="red",lwd=2)

res_vb_males<-residuals(vbgfmales)
fit_vb_males<-fitted(vbgfmales)
plot(res_vb_males~fit_vb_males,pch=19,col="seagreen",xlab="Fitted Values",ylab="Residual Values",main="Hake Male Growth Residuals Plot")

res_vb_females<-residuals(vbgffemales)
fit_vb_females<-fitted(vbgffemales)
plot(res_vb_females~fit_vb_females,pch=19,col="seagreen",xlab="Fitted Values",ylab="Residual Values",main="Hake Female Growth Residuals Plot")

Maturity

A maturity data subset was created using the age and gonad stage factors. A maturity function is then used to make a common fisheries model: an asymptotic curve. This model is then optimized using the NLS function and applyed to the plot. a50 = age at first capture (age at 50%), a95 = age between 50% - 95%. It is once again necessary to use the unique ages to increase accuracy. The plot then displays the prportion of hake that reach maturity at a certain age. Residual error is not required when plotting maturity due to the use of the function.

datMat<-subset(hake,!is.na(hake$age)&!is.na(hake$gonad_stage))
mature<-function(x){
  tot<-length(x)
  mat<-length(x[x>1])
  prop<<-mat/tot 
}
n<-max(datMat$age)
mat<-data.frame(age=1:n,prop=NA)
for(i in 1:n){
  dum<-subset(datMat,age==i&!is.na(gonad_stage))
  mat$prop[i]<-ifelse(nrow(dum)>0,mature(dum$gonad_stage),NA)
}
fitMat<-nls(prop~1/(1+19^((afif-age)/ani)),data=mat,start=c(afif=2,ani=0.5))
plot(prop~age,data=mat,main="Maturity Hake Combined",xlab="Age (Years)",ylab="Proportion of Hake",pch=19,col="green3")
uniqueages_mat<-seq(range(datMat$age)[1],range(datMat$age)[2],0.1)
text(x=10, y=0.5, paste("a50=", round(coef(fitMat)[1],3),sep=""))
text(x=15, y=0.5, paste("a95=", round(coef(fitMat)[2],3),sep=""))
lines(predict(fitMat,data.frame(age=uniqueages_mat))~uniqueages_mat,col="red3",lwd=2)

Males and Females The same maturity function and asymptotic curve are applied to both males & females. We see the consistent trend of females have longer lifespans compared to males.

malesMat<-subset(datMat,sex =="1")
nM<-max(malesMat$age)
matM<-data.frame(age=1:nM,prop=NA)
for(i in 1:nM){
  dumM<-subset(malesMat,age==i&!is.na(gonad_stage))
  matM$prop[i]<-ifelse(nrow(dumM)>0,mature(dumM$gonad_stage),NA)
}
fitMat_males<-nls(prop~1/(1+19^((afif-age)/ani)),data=matM,start=c(afif=2,ani=0.5))
plot(prop~age,data=matM,main="Maturity Hake (Males)",xlab="Age (Years)",ylab="Proportion of Hake",pch=19,col="green3")
text(x=7, y=0.5, paste("a50=", round(coef(fitMat_males)[1],3),sep=""))
text(x=11, y=0.5, paste("a95=", round(coef(fitMat_males)[2],3),sep=""))
uniqueages_mat_males<-seq(range(malesMat$age)[1],range(malesMat$age)[2],0.1)
lines(predict(fitMat_males,data.frame(age=uniqueages_mat_males))~uniqueages_mat_males,col="red3",lwd=2)

femalesMat<-subset(datMat,sex =="2")
nF<-max(femalesMat$age)
matF<-data.frame(age=1:nF,prop=NA)
for(i in 1:nF){
  dumF<-subset(femalesMat,age==i&!is.na(gonad_stage))
  matF$prop[i]<-ifelse(nrow(dumF)>0,mature(dumF$gonad_stage),NA)
}
fitMat_females<-nls(prop~1/(1+19^((afif-age)/ani)),data=matF,start=c(afif=2,ani=0.5))
plot(prop~age,data=matF,main="Maturity Hake (Females)",xlab="Age (Years)",ylab="Proportion of Hake",pch=19,col="green3")
text(x=10, y=0.5, paste("a50=", round(coef(fitMat_females)[1],3),sep=""))
text(x=15, y=0.5, paste("a95=", round(coef(fitMat_females)[2],3),sep=""))
uniqueages_mat_females<-seq(range(femalesMat$age)[1],range(femalesMat$age)[2],0.1)
lines(predict(fitMat_females,data.frame(age=uniqueages_mat_females))~uniqueages_mat_females,col="red3",lwd=2)

Mortality

For mortality, the histogram shows the trend of the hake population mortality frequency at each age. Peak mortality is observed to occur between ages 5-7. We then log-transform the mortality data. We see it is successful so far, as the trend shape is consistent across the histogram and dot plot.

datMort<-subset(hake,!is.na(hake$age))
dt<-hist(datMort$age,main="Hake Mortality",xlab="Age (Years)",ylab="Frequency",col="goldenrod",breaks=seq(from=min(datMort$age)-0.5,
                            to=max(datMort$age)+0.5,by=1))

ds<-data.frame(age=dt$mids,logfreq=log(dt$counts))
plot(logfreq~age,data=ds,main="Hake Mortality (log Frequency)",xlab="Age (Years)",ylab="Log (Frequency)",pch=19,col="goldenrod")

The log linear transformation assumes constant exploitation. To allow a straight trend line to be added, the graph must have more selectivity of age, by eliminating younger ages up to the highest age to the lowest age. After the linear model is added, we see a clear mortality linear model. Z = total mortality, M = natural mortality, F = fished mortality (Z-M).

d_final<-ds[5:20,]
plot(logfreq~age,data=d_final,xlab="Age",ylab="Log (Frequency)",main="Hake Mortality (log Frequency)",pch=19,col="goldenrod")
text(x=6, y=3, paste("Z = 0.23",sep=""))
text(x=7, y=2.5, paste("M = 0.19",sep=""))
text(x=8, y=2, paste("F = 0.002",sep=""))
resMort<-lm(logfreq~age,data=d_final)
lines(d_final$age,predict(resMort),col="red")

YPR Model

The yield per recruit (YPR) model will be used to identify, explore and evaluate the current status of the hake stock.

Using the previous calculations, the resulting parameters will now be used for the YPR model.

param=list()
# No. recruits (age0)
param$rec=1000
#Growth
param$Linf=113.65
param$K=0.163
param$t0=-1.209
#Weight at Length
param$a=0.002
param$b=3.404
#Maturity
param$af=4.233
param$ato=1.933
#Mortality
param$M=0.19

Now the YPR model is developed using a function that includes fishing mortality, age at first capture, a95, and the parameters selected.

YPR_Model<-function(Fm,afc,sel2,param){
  with(param,{
    Age=1:26
    LengthAtAge<<-Linf*(1-exp(-K*(Age-t0)))
    WeightAtLength<<-a*LengthAtAge^b
    PropMature<<-1/(1+19^((af-Age)/ato))
    SelectAtAge<<-1/(1+19^((afc-Age)/sel2))
    
    NAtAge<-rep(0,length(Age))
    for(i in 2:length(Age)){
      NAtAge[1]<-rec*exp(-(M+Fm*SelectAtAge[1]))
      NAtAge[i]<-NAtAge[i-1]*exp(-(M+Fm*SelectAtAge[i]))
    }
    NumberAtAge<<-NAtAge
    
    NDying<-rep(0,length(Age))
  for(i in 2:length(Age)){
    NDying[1]<-rec*(1-exp(-(M+Fm*SelectAtAge[1])))
    NDying[i]<-NAtAge[i-1]*(1-exp(-(M+Fm*SelectAtAge[i])))
  }
  NumberDying<<-NDying
  
  NumberCaught<<-(Fm*SelectAtAge)/(M+Fm*SelectAtAge)*NumberDying
  CatchAtAge<<-NumberCaught*WeightAtLength
  YieldPerRecruit<-sum(CatchAtAge)/rec
  return(YieldPerRecruit)
  })
}

Now, we use the model with our estimated fishing mortality to predict the yield per recruit, assuming the initial age of first capture is our calculated figure of 4.32. The YPR result of 318 suggests the stock is in a positive trend of reproduction and the hake stock is not under serious fishing pressure.

YPR_Model(Fm=0.002,afc=4.23,sel2=1.93,param)
## [1] 45.36846

Optimization

Firstly, we optimize the Fmax, or the highest rate of fishing mortality that maximizes the yield from the hake fishery. This increases to from 0.002 to 0.131.

Next, we optimize the age at first capture to find the highest age that maximizes the YPR, using the new Fmax. This is discovered to be 9.7 years, increasing from 4.23 years.

optimize(YPR_Model,afc=4.23,c(0,2),sel2=1.93, param=param,maximum = TRUE)
## $maximum
## [1] 0.4388486
## 
## $objective
## [1] 1782.496
optimize(YPR_Model,Fm=0.131,c(1,26),sel2=1.93, param=param,maximum = TRUE)
## $maximum
## [1] 3.912454
## 
## $objective
## [1] 1426.823

Role of AFC

Now we plot the effect of the 2 differing age at first captures calculated: the original AFC at maturity, and the optimized AFC.

fishing<-seq(0,2,0.1)
YPR<-1
for(i in 1:length(fishing)){
  YPR[i]<-YPR_Model(Fm=fishing[i],afc=4.23,sel2=1.93,param)
}
plot(YPR~fishing,type="l",col="purple",lwd=2,ylab="Yield per Recruit",xlab="Fishing Mortality",main="Hake Fishing Effort",ylim=c(0,3000))
YPR_2<-1
for(i in 1:length(fishing)){
  YPR_2[i]<-YPR_Model(Fm=fishing[i],afc=9.70,sel2=1.93,param)
}
points(YPR_2~fishing,type="l",lty=2,lwd=2,col="red")
point1 <- YPR_Model(Fm = 0.002, afc = 9.7, sel2 = 1.93, param)
point2 <- YPR_Model(Fm = 0.131, afc = 9.7, sel2 = 1.93, param)
points(0.002, point1, pch=19, col="69")   # fished mortality
points(0.131, point2, pch=17, col="green") # maximum YPR

# Updated legend
legend("topright",
       legend = c("AFC = 4.23", "Optimal AFC = 9.70", 
                  "Fm = 0.002", "Fmax = 0.131"),
       col = c("purple", "red", "69", "green"),
       lty = c(1, 2, NA, NA),
       pch = c(NA, NA, 19, 17),
       lwd = c(2, 2, NA, NA))

Spawner Stock Biomass In order to consider sustainability, we need to know how many fish are in the water, not just how many we have caught. This is based on the spawner stock biomass (SSB). Firstly we develop a function to create the SSB model. This is our Spawner per Recruit (SPR) model.

SPR_Model <- function(Fm, afc, sel2, param){
  with(param, {
    
    Age=1:26
    
    LengthAtAge <<- Linf*(1-exp(-K*(Age-t0)))
    WeightAtLength <<- a*LengthAtAge^b
    PropMature <<- 1/(1+19^((af - Age)/ato)) 
    SelectAtAge <<- 1/(1+19^((afc- Age)/sel2))
    
    NAtAge <- rep(0, length(Age))  
    for(i in 2:length(Age)){
      NAtAge[1] <- rec*exp(-(M+Fm*SelectAtAge[1]))
      NAtAge[i] <- NAtAge[i-1]*exp(-(M+Fm*SelectAtAge[i]))
    }
    NumberAtAge <<- NAtAge
    
    NDying <- rep(0, length(Age))
    for(i in 2:length(Age)){
      NDying[1] <- rec*(1-exp(-(M+Fm*SelectAtAge[1])))
      NDying[i] <- NAtAge[i-1]*(1-exp(-(M+Fm*SelectAtAge[i])))
    }
    NumberDying <-- NDying
    
    NumberCaught <<- (Fm*SelectAtAge)/(M+Fm*SelectAtAge)*NumberDying 
    
    CatchAtAge <<- NumberCaught*WeightAtLength 
    
    BiomassAtAge <<- NumberAtAge*WeightAtLength
    
    SSBAtAge <<- NumberAtAge*WeightAtLength*PropMature
    
    SpawnerPerRecruit <- sum(SSBAtAge)/rec
    
    return(SpawnerPerRecruit)
    
  })
  
}

VirginSSB: We then create a baseline spawner stock biomass when there is no fishing pressure, assuming the optimum AFC of 9.7.

We then find the minimum sustainable spawner stock biomass, which is 40% of the virginSSB. This is the accepted fisheries benchmark for sustainable SPR, to allow for 40% of the original reproductive capacity to remain untouched.

We then caclulate the current SSB to see whether a fishing pressure of 0.002 at an AFC of 9.7 is sustainable. It is less than the virgin biomass, leaving 98.5% of the virgin biomass untouched. This is very sustainable.

VirginSSB<-SPR_Model(Fm=0,afc=9.7,sel2=1.93,param=param)
VirginSSB
## [1] 20885.99
SustainableSSB<-0.4*VirginSSB
SustainableSSB
## [1] 8354.395
SSBnow<-SPR_Model(Fm=0.002,afc=9.7,sel2=1.93,param=param)
SSBnow
## [1] 20754.53
SSBnow/VirginSSB
## [1] 0.9937059

Optimizing SPR40: This will allow us to calculate the maximum Sustainable Yield (MSY) for the hake fishery. First, we calculate our SPR40 function for the model.

SPR_40 <- function(Fm, afc, sel2, param, SustainableSSB){
  
  with(param, {
    
    Age=1:26
    
    LengthAtAge <<- Linf*(1-exp(-K*(Age-t0)))
    WeightAtLength <<- a*LengthAtAge^b
    PropMature <<- 1/(1+19^((af - Age)/ato)) 
    SelectAtAge <<- 1/(1+19^((afc- Age)/sel2)) 
    
    NAtAge <- rep(0, length(Age)) 
    for(i in 2:length(Age)){
      NAtAge[1] <- rec*exp(-(M+Fm*SelectAtAge[1]))
      NAtAge[i] <- NAtAge[i-1]*exp(-(M+Fm*SelectAtAge[i]))
    }
    NumberAtAge <<- NAtAge
    
    NDying <- rep(0, length(Age))
    for(i in 2:length(Age)){
      NDying[1] <- rec*(1-exp(-(M+Fm*SelectAtAge[1])))
      NDying[i] <- NAtAge[i-1]*(1-exp(-(M+Fm*SelectAtAge[i])))
    }
    NumberDying <-- NDying
    
    NumberCaught <<- (Fm*SelectAtAge)/(M+Fm*SelectAtAge)*NumberDying 
    
    CatchAtAge <<- NumberCaught*WeightAtLength 
    
    BiomassAtAge <<- NumberAtAge*WeightAtLength
    
    SSBAtAge <<- NumberAtAge*WeightAtLength*PropMature
    
    SpawnerPerRecruit <- sum(SSBAtAge)/rec
    
    SustainableSSB = SustainableSSB
    
    Difference<<-(SpawnerPerRecruit-SustainableSSB)^2
    
    return(Difference)
    
  })
  
}

Optimized MSY The returns using the SPR40 model, assuming current fishing pressure and optimum AFC. Now we optimize the fishing pressure for MSY, ideally with a minimal difference between this and the sustainable SSB.

The difference is very low between the SSB and the new fishing pressure (0.0813264). Our new fishing pressure/mortality rate to achieve FMSY is 0.1875756. When this is applied to our original SPR model, we finally end up with our MSY of 63893.02. This is very very similar to our spawner stock biomass of 63892.74, separated by 0.28.

SPR_40(Fm=0.002,afc=9.7,sel2 = 1.93, param=param, SustainableSSB = SustainableSSB)
## [1] 153763293
optimize(SPR_40,c(0,2),afc=9.7,sel2=1.93,param=param, SustainableSSB)
## $minimum
## [1] 1.742924
## 
## $objective
## [1] 0.000183843
SPR_Model(Fm=0.1875756,afc=9.7,sel2=1.93,param=param)
## [1] 14183.94
SustainableSSB
## [1] 8354.395

Climate Modelling

To assess how rising sea temperatures might affect the future productivity of the New Zealand hake fishery, a modified yield-per-recruit (YPR) model was developed incorporating climate-related biological shifts. Based on predictions from multispecies size-spectrum models (Datta et al., 2024), key parameters were adjusted to reflect likely responses of hake to warming: increased growth rate (K), reduced maximum length (L∞), earlier maturity (af), and higher natural mortality (M). Under this “warming scenario,” the YPR curve demonstrated a noticeable reduction in yield at all fishing mortalities, with a lower Fmax and maximum yield. These results suggest that climate change could diminish the productivity of the hake fishery, reducing sustainable catch levels unless fishing effort is adjusted accordingly.

# Simulate warming scenario
param_warm <- param
param_warm$K <- 0.18      # Slightly faster growth
param_warm$Linf <- 105    # Lower asymptotic length
param_warm$af <- 3.5      # Earlier age at maturity
param_warm$ato <- 1.5     # Quicker transition to maturity
param_warm$M <- 0.22      # Higher natural mortality

# Plot YPR curves for both scenarios
fishing <- seq(0, 2, 0.05)

# Baseline
YPR_baseline <- sapply(fishing, function(F) YPR_Model(Fm = F, afc = 9.7, sel2 = 1.93, param = param))

# Warming scenario
YPR_warm <- sapply(fishing, function(F) YPR_Model(Fm = F, afc = 9.7, sel2 = 1.93, param = param_warm))

# Plot
plot(fishing, YPR_baseline, type = "l", lwd = 2, col = "purple", ylab = "Yield per Recruit", ylim=c(0,3000),
     xlab = "Fishing Mortality (F)", main = "YPR Comparison: Baseline vs. Warming Scenario")
lines(fishing, YPR_warm, lwd = 2, col = "red", lty = 2)
legend("topright", legend = c("Baseline", "Warming scenario"),
       col = c("purple", "red"), lty = c(1, 2), lwd = 2)

Now we use the new warming scenario model to figure our our new parameters. The warming virgin biomass of 13775,a nd the Sustainable SSB under the warming scenarios

VirginSSB_warm <- SPR_Model(Fm = 0, afc = 9.7, sel2 = 1.93, param = param_warm)
SustainableSSB_warm <- 0.4 * VirginSSB_warm
VirginSSB_warm
## [1] 13775.13
SustainableSSB_warm
## [1] 5510.053

The optimized fishing mortality under warming. This drastically increases from 0.002 to 1.99

FMSY_warm <- optimize(SPR_40,
                      c(0, 2),
                      afc = 9.7,
                      sel2 = 1.93,
                      param = param_warm,
                      SustainableSSB = SustainableSSB_warm)$minimum
FMSY_warm
## [1] 1.99994

The Yield of the warming MSY is 6829, which is 10% of the baseline MSY.

SPR_Model(Fm = 1.994, afc = 9.7, sel2 = 1.93, param = param_warm)
## [1] 6829.313
Discussion

This stock assessment of Merluccius australis (southern hake) in New Zealand reveals a population currently under minimal fishing pressure, with strong indications of sustainability. At a current estimated fishing mortality of F = 0.002, and an age at first capture (AFC) of 9.7 years, the yield per recruit (YPR) is 318, and 98.5% of the virgin spawning biomass remains intact. This level of spawning potential is well above the accepted 40% sustainability threshold (SPR40), indicating that the fishery is operating conservatively. These model-based findings align closely with the 2024 Fisheries NZ Plenary Report, which states that all three hake stocks (HAK 1, HAK 4, and HAK 7) are currently above their respective management targets of B₄₀.

Optimizing fishing pressure and selectivity further strengthened this conclusion. An optimal fishing mortality (Fmax = 0.131) and delayed AFC of 9.7 years produced a peak YPR and a maximum sustainable yield (MSY) estimate of ~63,893, nearly matching the calculated spawning stock biomass (SSB) under that regime. This near-perfect yield-to-biomass alignment suggests that the current fishery has considerable headroom before approaching biological overfishing limits, consistent with previous observations that hake in New Zealand are typically caught as bycatch in the hoki and blue whiting fisheries rather than as a primary target (Horn, 2015; Arancibia et al., 2015).

However, this robust current status must be viewed in the context of a rapidly changing marine environment. Hake are a deepwater species whose productivity is strongly tied to environmental stability and prey availability in regions like the Chatham Rise—a biologically productive but climatically sensitive area. According to Pinkerton et al. (2021) and Datta et al. (2024), deep-sea ecosystems in New Zealand, including the Sub-Antarctic and Chatham Rise, are likely to experience declines in biomass and fishery yield under projected warming scenarios, with effects intensified by fishing pressure.

To explore these dynamics, this study incorporated a warming scenario into the YPR and SPR models. Based on expected physiological responses to ocean warming (e.g. increased growth rate, reduced maximum length, earlier maturity, and higher natural mortality), the altered parameters led to dramatic declines in productivity. The MSY under the warming scenario dropped to 6,829, representing a 90% decline compared to baseline, and only achievable under extremely high fishing mortality (FMSY = 1.99), which would be biologically unrealistic and likely unsustainable.

These results are consistent with recent multispecies ecosystem modelling of New Zealand waters using mizer and therMizer, which show that deepwater species are particularly vulnerable to warming, especially under continued fishing pressure (Datta et al., 2024). Furthermore, observed shifts in species distributions and depth ranges in New Zealand (e.g., Sutton & Bowen, 2014; Pinkerton et al., 2021) support concerns that hake may become less accessible or productive in current fishing grounds.

Management actions should therefore be forward-looking. The benefits of delaying AFC to 9.7 years were evident even under warming, suggesting that gear modifications or regulatory minimum size limits could be useful tools for future-proofing the fishery. Maintaining high spawning biomass will buffer the population against recruitment failure or stochastic environmental events. Additionally, adaptive harvest strategies—such as tiered TACs based on environmental indices, or real-time catch monitoring—may be necessary as biological parameters shift.

In summary, the New Zealand hake fishery is presently in a healthy state, with ample reproductive capacity and conservative exploitation. However, this resilience is not guaranteed in the face of warming seas. This study shows that the productivity and sustainability of the fishery could deteriorate rapidly under even modest climate-induced biological shifts. By embedding climate sensitivity into stock assessments and continuing to monitor hake biology and distribution, fisheries managers can maintain this success story well into the future.

References

Anderson, O. F., Edwards, C. T., & Ballara, S. L. (2019). Non-target fish and invertebrate catch and discards in New Zealand hoki, hake, ling, silver warehou, and white warehou trawl fisheries from 1990–91 to 2016–17 (p. 117). Ministry for Primary Industries, Manatū Ahu Matua.

Arancibia, H., et al. (2015). An overview of hake and hoki fisheries: analysis of biological, fishery, and economic indicators.

Datta, S., Beran, H., & Rogers, A. (2024). The impacts of warming on shallow and deep-water fisheries in New Zealand. Preprint, Earth’s Future.

Fisheries New Zealand. (2024). Fisheries Assessment Plenary May 2024: Stock Assessments and Yield Estimates.

Horn, P. L. (2015). Southern hake (Merluccius australis) in New Zealand: Biology, fisheries and stock assessment. In Hakes: Biology and exploitation (pp. 101–125).

Pinkerton, M. H., Forman, J. S., & Tait, A. (2021). Climate Change and New Zealand Fisheries: Impacts, Risks and Management Options. NIWA Science and Technology Series.

Sutton, P. J. H., & Bowen, M. (2014). Ocean temperature change around New Zealand over the last 36 years. New Zealand Journal of Marine and Freshwater Research, 48(2), 273–284.