- Email: [email protected]

Contents lists available at ScienceDirect

Fisheries Research journal homepage: www.elsevier.com/locate/fishres

Improving management advice through spatially explicit models and sharing information Carlos Montenegro a,∗ , Mark N. Maunder b , Maximiliano Zilleruelo a a b

Instituto de Fomento Pesquero, Evaluacion de Recursos, Blanco 839, Valparaíso, Chile Inter-American Tropical Tuna Comisión, 8604 La Jolla Shores Drive, La Jolla, CA 92037-1508, USA

a r t i c l e

i n f o

Article history: Received 23 April 2009 Received in revised form 23 July 2009 Accepted 24 July 2009 Keywords: Fisheries management Heterocarpus reedi Shrimp Spatial model Stock assessment

a b s t r a c t Recent assessments of Chilean shrimp, Heterocarpus reedi, in central Chile have been conducted separately for the northern and southern zones of the ﬁshery and treating them as two separate stocks. However, it is not clear whether H. reedi of the two zones interact with one another or whether they share similar characteristics. Such knowledge is necessary to determine whether they should be modeled as separate “stocks” or as a single stock. This has motivated the use of the Pella–Tomlinson model to test whether there are spatial differences in the population dynamics of H. reedi in the two zones and whether sharing information between the zones improves management advice. We test if it is better, from a stock assessment point of view, to model the stock as one unit in the whole area, or as two separate stocks. In the single-stock model, we sum the catch data of both zones, but each catch-per-unit-of-effort index is ﬁt as a separate data set, using a joint likelihood. Under the single-stock hypothesis, the best model ﬁt was the symmetric production function (i.e. the Schaefer model for which the biomass that supports maximum sustainable yield as a proportion of carrying capacity (BMSY /B0 ) = 0.5), with different catchability coefﬁcients for each CPUE index, but a shared standard deviation of the log-normal likelihood function. Under the two-stock hypotheses, both catch and CPUE data were separated for each zone in the model. In this case, the best model ﬁt is also the one with symmetrical production curve, and the only parameter that differed between the zones was B0 . However, B0 per unit of habitat was similar for the two zones. Also, the precision of estimated management quantities was improved by modeling the appropriate spatial structure and sharing information among zones. The results suggest that the demographic parameters are similar for the two zones. It appears that the main difference between the two zones is the exploitation history, with the catch in the southern zone being reduced earlier than in the northern zone and consequently the biomass in the southern zone increased earlier than in the northern zone. This implies that local depletion can occur in this stock and that differences in management among zones may require explicitly modeling sub-stocks in the assessment of this and other species. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The Chilean shrimp (Heterocarpus reedi, Decapoda, Pandalidae) has a wide distribution off the coast of Chile. It inhabits the bottom of the continental shelf and upper slope from about 21◦ 30 S to about 38◦ 43 S, at depths between 150 and 600 m. Commercial exploitation started during the 1960s, making it one of the oldest crustacean ﬁsheries in Chile. The ﬁshery extends between region II and region VIII (Region de Atacama to Region del Bío Bío), but the main ﬁshing grounds are located between 27◦ S and 37◦ S (Region III to Region VIII). The ﬁshery administration (Subsecretariat of Fishery) declared the ﬁshery to be fully exploited in 1995,

∗ Corresponding author. Tel.: +56 32 2151404; fax: +56 32 2151645. E-mail address: [email protected] (C. Montenegro). 0165-7836/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.ﬁshres.2009.07.006

and established a total allowable catch quota of 10,000 t for 1996. Spatially, the management of the resource is based on two main zones: the northern zone, from latitude 26◦ 03 S to 32◦ 12 S (including Regions II to IV), and the southern zone, from latitude 32◦ 12 S to 38◦ 28 S (Regions V to VIII, Fig. 1). As landings showed a precipitous decrease during 1997–2000, individual quotas were introduced in 2001 and the total catch quota was split in a temporal and spatial fractioning that included the complete closure of the southern zone. Since then, alternate closures to regions in the southern zone have been applied annually, while a total catch quota has been maintained at around 5000 t/year. The landings have been around 80% of total catch quota with a total of 4215 t in 2006 and 4340 t in 2007. No spatial restrictions have been applied since 2006. The management of the ﬁshery includes annual stock assessments to estimate the status of the resource, analyze the implications of different management actions, and determine the total allowable catch.

192

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

Fig. 2. Yield (surplus production) as a proportion of MSY at different levels of biomass (as a fraction of virgin biomass) for different values of shape parameter.

mine if the virgin biomass (B0 ), the biomass at the beginning of the modeling period (1989) as a proportion of B0 (p) or the production rate at maximum production (r), differs between the two zones that H. reedi inhabits. Furthermore, we test whether the data suggest a better ﬁt to a symmetric production curve (where the biomass at maximum sustainable yield (MSY) as a ratio of the unexploited biomass (BMSY /B0 = 0.5) or an asymmetric (BMSY /B0 = / 0.5) one). We use the Gilbert’s (1992) formulation of the Pella–Tomlinson model (see Maunder and Starr, 1995). The general formulation for a surplus-production model is (in difference equation or discrete form): Bt+1 = Bt + g(Bt ) − Ct

Fig. 1. Geographic distribution of the ﬁshery of Heterocarpus reedi off the coast of Chile. Fishery zones are also indicated.

These assessments have used an age-structured catch-at-length model. However, it is not clear whether H. reedi in the two zones interact or if they share similar characteristics. Such knowledge is necessary to determine whether they should be modeled as separate “stocks” or as a single stock. The goal of this paper is to determine whether there are spatial differences in the population dynamics of H. reedi of the coast of Chile. Therefore, we test whether it is better to model the stock as one unit in the whole area, or as two separate stocks. Speciﬁcally, we investigate whether modeling the spatial structure and, following Polovina (1989), if sharing information between zones improves management advice. 2. Materials and methods The Pella–Tomlinson (1969) model is used to represent the dynamics of the stocks. We use “model” to deﬁne the shape parameter assumption of the Pella–Tomlinson model, “hypothesis” to deﬁne one or two stocks, and “scenarios” to deﬁne which parameters are shared between the two stocks. The data used to evaluate stock-structure hypotheses, models, and scenarios are the series of landings and the standardized catch-per-unit-of-effort (CPUE) from 1989 to 2007. The landings from 1989 to 2007 were obtained from the ofﬁcial statistics of Servicio Nacional de Pesca of Chile (www.sernapesca.cl), while the standardized CPUE was estimated using a GLM (Generalized Linear Model, McCullagh and Nelder, 1989). We use three versions of a surplus-production model, deﬁned by the value of the Pella–Tomlinson model shape parameter, and we test several scenarios. The scenarios differed as to which parameters were shared by the two zones and, consequently, as to the number of parameters estimated. This analysis is used to deter-

(1)

where Bt+1 is the exploitable biomass at the beginning of year t + 1, g(Bt ) is the production of biomass as a function of biomass at the start of the year t, and Ct is the biomass caught during year t. The production function of biomass, in the formulation of Gilbert (1992) is: r g(Bt ) = (1/m) − 1

Btm B0m−1

− Bt

(2)

Eq. (2) is not valid when m = 1. The Pella–Tomlinson model permits consideration of the Schaefer (1954) and Fox (1970) models as special cases of the generalized form. The equation has three parameters: B0 , the virgin biomass, m, the shape parameter of the production curve, and r, the production rate at maximum production. The ratio BMSY /B0 is calculated using the expression: BMSY 1 = B0 m1/(m−1)

(3)

The production parameter differs from the typical deﬁnition of the rate when the biomass approaches zero. r=

g(BMSY ) BMSY

(4)

We test three alternative models, based in the shape parameter: Model

BMSY /B0

m

Model 1 Model 2 Model 3

0.5 0.3 –

2 0.68 Free parameter

Model 1 corresponds to a classical Schaefer model with a symmetrical production curve (see Fig. 2). Model 2 ﬁxes the shape parameter based on the biology (natural mortality, growth rate, and steepness of the stock–recruitment relationship) of H. reedi, as suggested by Maunder (2003). Roa and Ernst (1996) report the growth rate for H. reedi to be 0.174 (year−1 ) for

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

193

Table 1 Scenarios used for the single-stock hypothesis and AIC values for the three models. Scenario

Number of parameters

Parameters

AIC BMSY /B0

Base model Scenario 1 Scenario 2 Scenario 3

5 6 6 7

B0 , r, p, 1 , q1 B0 , r, p, 1 , 2 , q1 B0 , r, p, 1 , q1 , q2 B0 , r, p, 1 , 2 , q1 , q2

Model 1 0.5

Model 2 0.3

Model 3 Estimated

−27.81 −26.98 −28.57 −27.49

−27.15 −25.18 −27.87 −25.88

−26.78 −25.11 −27.53 −25.74

Bold values means best ﬁts.

females and 0.199 (year−1 ) for males. Canales et al. (1999; http://www.ﬁp.cl/pdf/informes/infﬁnal%2097-24.pdf) report natural mortality of 0.28 (year−1 ). Currently, there is no available information about the stock–recruitment relationship for this resource. Therefore, based on Table 1 of Maunder (2003) (with natural mortality = 0.3 (year−1 ), growth rate = 0.2 (year−1 ), and a steepness of the stock–recruitment relationship = 0.75), we decided to ﬁx the shape parameter at 0.68 for model 2, which corresponds to a ratio BMSY /B0 of 0.3. Model 3 allows the shape parameter to be estimated along with the rest of the model parameters. The single-stock hypothesis is based on a model that represents the dynamics of the stock as one unit, for which the biomass at time t + 1 is a function of the biomass at time t, the parameters of the production function, and the total landings. The exploitable biomass at time t + 1 is: Bt+1 = Bt +

r (1/m) − 1

Btm B0m−1

− Bt

2 −

Ct,s

(5)

s=1

where Ct,s is the landings per year (t) and zone (s). The biomass at the beginning of the ﬁrst year is a proportion of the virgin biomass, so we have an extra parameter (p): B1989 = pB0

(6)

The observation model corresponding to the CPUE data is: It,s = qs Bt evt,s

(7)

where It,s is the observed standardized CPUE, which was estimated from the application of a GLM procedure, qs is the zone-speciﬁc catchability coefﬁcient, which was also estimated in the ﬁtting procedure, and t,s is a normally distributed observation error (∼N(0, s2 )).

The model was ﬁtted by maximum likelihood, assuming a log-normal distribution for the standardized CPUE. An index of standardized CPUE is available for each zone, and both indices are included in the likelihood function (negative log-likelihood function, omitting constants), as follows: ln(L) =

T 2

ln(ˆ s ) +

(ln(It,s ) − ln(ˆqs Bˆ t ))

2

(8)

2ˆ s2

s=1 t=1

The standard deviation of the likelihood function (ˆ s ) is also estimated in the ﬁtting procedure. We used a forward stepwise algorithm based on the AIC criteria (Hilborn and Mangel, 1997) to determine which model parameters should be shared between the two zones. For the single-stock hypothesis, the CPUE parameters (q and ) are considered to differ for the two zones. This required running the three models previously described under four scenarios (Table 1), where q1 and q2 are the catchability parameters for the CPUE in the north and south, respectively, and 1 and 2 are the standard deviation parameters for the CPUE in the north and south, respectively. In the case of Model 3, we also estimate the shape parameter of the production curve (m) increasing the number of parameters by one (Table 1). The two-stock hypothesis is based on a model that represents the dynamics of the stock for each zone. The biomass at time t, for stock s is: Bt+1,s

rs = Bt,s + (1/ms ) − 1

ms Bt,s ms −1 B0,s

− Bt

− Ct,s

(9)

Again, the biomass at the beginning of the ﬁrst year (for each stock), is a proportion of the virgin biomass: B1989,s = ps B0,s

(10)

Table 2 Scenarios used for the two-stock hypothesis and AIC values for the three models. Scenario

Number of parameters

Parameters

AIC BMSY /B0 Model 1 0.5

Base model Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5 Scenario 6 Scenario 7 Scenario 8 Scenario 9 Scenario 10 Scenario 11 Bold values means best ﬁts.

5 6 6 6 6 6 7 7 7 7 8 10

B01 , r1 , p1 , 1 , q1 B01 , B02 , r1 , p1 , 1 , q1 B01 , r1 , r2 , p, 1 , q1 B01 , r1 , p1 , p2 , 1 , q1 B01 , r1 , p1 , 1 , 2 , q1 B01 , r1 , p1 , 1 , q1 , q2 B01 , B02 , r1 , r2 , p, 1 , q1 B01 , B02 , r1 , p1 , p2 , 1 , q1 B01 , B02 , r1 , p1 , 1 , 2 , q1 B01 , B02 , r1 , p1 , 1 , q1 , q2 B01 , B02 , r1 , p1 , 1 , q1 , q2 , m2 B01 , B02 , r1 , r2 , p1 , p2 , 1 , s2 , q1 , q2

−19.98 −46.58 −41.52 −44.44 −20.58 −35.43 −44.74 −44.61 −44.58 −44.88 – −39.13

Model 2 0.3 −22.04 −43.78 −39.02 −43.21 −23.03 −35.66 −42.08 −42.07 −41.90 −42.07 – –

Model 3 Estimated −21.26 −45.65 −40.30 −43.27 −22.54 −34.08 −43.68 −44.13 −43.78 −43.77 −44.32 –

194

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

The observation model corresponding to the CPUE data in each zone is: It,s = qs Bt,s evt,s

(11)

The likelihood function (negative log-likelihood function, omitting constants) is the same as (8) except the biomass is stock-speciﬁc. For the two-stock hypothesis, the forward stepwise algorithm required running the three models previously described under 10 scenarios, differing in whether the population dynamics model parameters or the CPUE index parameters differed between the two stocks (where 1 indicates the northern stock and 2 indicates the southern stock). In the case of Model 3, we also estimate the shape parameter of the production curve (m), which is assumed the same for each stock, increasing the number of parameters by one (Table 2). For the comparison of models and scenarios, we use the Akaike information criterion (AIC), which is deﬁned as two times the negative log-likelihood, plus two times the number of parameters (Hilborn and Mangel, 1997).

AIC = −2 ln L /data

+ 2

(12)

AD Model Builder (developed by David Fournier and now freely available from admb-project.org) is used to implement the model. The parameters of the model were estimated in log scale (to avoid negative values), so the variance of the model’s parameters and the variance of the derived quantities (e.g. biomass and MSY) were estimated using the delta method (Seber, 1973). 3. Results We test the two spatial structure hypothesis proposed in Section 2 (single stock versus two stocks) and three formulations for the production function (m = 2, m = 0.68 and m = free parameter). For each of these hypotheses and models, we ﬁt scenarios of ascending complexity. For the case of the models under the hypothesis of a single stock, the best model ﬁt was Model 1 (production function with m = 2, BMSY /B0 = 0.5), with different catchability coefﬁcients for each CPUE index. For this model, we estimate virgin biomass of 63,556 t, the biomass at the beginning of 1989 of 16,562 t (26% of B0 ), and a rate of production at maximum production (r) of 0.31. For this scenario, the biomass at maximum sustainable yield (BMSY ) is 31,778 t, with MSY of 9857 t (Table 3). This model has an AIC = −28.57, −ln(L) = −40.57 and six parameters (B0 , p, r, , q1 , q2 ). Model 2 (with m = 0.68; BMSY /B0 = 0.3) estimated virgin biomass several orders of magnitude greater than that for Model 1, with a low rate of production (r < 0.12). In the case of Model 3 (m = free parameter) the best model ﬁt (AIC = −27.53) estimates a virgin biomass of 39,145 t, the biomass at the beginning of 1989 of 14,218 t and a rate of production at maximum production (r) of 0.4. In Model

Table 3 Parameters, management quantities, standard errors, coefﬁcient of variation (CV%) and 95% conﬁdence limits for the best models of the single- and two-stock hypotheses. Parameter

Lower bound

Upper bound

Best model of the single-stock hypothesis 63,556 28,418 45% B0 p 0.26 0.07 28% r 0.31 0.09 30% 0.20 0.02 12% 7.48E−06 1.61E−06 22% q1 q2 8.79E−06 1.85E−06 21% BMSY 31,778 14,209 45% MSY 9857 4407 45%

Value

Std. Err

CV%

7857 0.119 0.126 0.151 4.32E−06 5.16E−06 3928 1219

119,255 0.402 0.494 0.242 1.06E−05 1.24E−05 59,628 18,495

Best model of the two-stock hypothesis 28,578 5968 B01 35,227 7491 B02 p 0.254 0.033 r 0.315 0.046 0.119 0.014 Q 1.63E−05 1.68E−06 14,289 2984 BMSY1 17,613 3745 BMSY2 MSY1 4500 325 MSY2 5547 440

16,880 20,545 0.190 0.225 0.092 1.30E−05 8440 10,273 3864 4685

40,276 49,909 0.319 0.405 0.147 1.96E−05 20,138 24,954 5136 6409

21% 21% 13% 15% 12% 10% 21% 21% 7% 8%

3, all scenarios suggest an asymmetrical production curve, with values of ratio BMSY /B0 between 0.55 and 0.65. For the case of the models under the hypothesis of two stocks, the best model ﬁt was again Model 1 (production function with m = 2, BMSY /B0 = 0.5), with a different virgin biomass for each zone. For this model, we estimate a total virgin biomass of 63,805 t (28,578 t for northern zone and 35,227 t for southern zone) and a total biomass at the beginning of 1989 of 16,213 t (7262 and 8951 t, each zone, 25% of B0 ); quantities very similar to the estimates of the best-ﬁt model of the single-stock hypothesis. Furthermore, the estimated r is 0.31 for both hypotheses. For this scenario, the biomass at MSY (BMSY ) is 14,289 t for northern zone and 17,613 t for the southern zone, with MSYs of 4500 and 5547 t, respectively (Table 3). This model has an AIC = −46.58, −ln(L) = −58.58 and six parameters (B01 , B02 , p, r, , q) and is the best model of both hypotheses, models, and scenarios. As for the single-stock hypothesis, Model 2 (with m = 0.68; BMSY /B0 = 0.3) estimated virgin biomass several orders of magnitude greater than Model 1 in several scenarios. In the case of Model 3 (m = free parameter) the best model ﬁt (AIC = −45.65) estimated a virgin biomass of 48,539 t (21,501 and 27,038 t each zone), the biomass at the beginning of 1989 of 14,218 t (6785 and 8532 t each zone) and an r of 0.36. As with the Model 3, almost all the scenarios suggest an asymmetrical production curve, with values of ratio BMSY /B0 between 0.4 and 0.65. The goodness of ﬁt of the models under the hypothesis of two stocks is better than the model under the hypothesis of one stock

Fig. 3. Observed CPUE (dots) and ﬁtted CPUE (line) in tons/trawling time from the best ﬁt of the single-stock hypothesis.

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

195

Fig. 4. Observed CPUE (dots) and ﬁtted CPUE (line) in tons/trawling time from the best ﬁt of the two-stock hypothesis.

Fig. 5. Results from the single-stock hypothesis. (a) Exploitable biomass estimated with the surplus-production model with the best ﬁt for the single-stock hypothesis (solid line, with 95% conﬁdence intervals), BMSY (horizontal solid line); (b) landings (solid line) and MSY (horizontal dashed line).

(AIC-based), which can be seen in the ﬁts to the observed CPUE (Figs. 3 and 4) and in the estimates of the standard deviations of the CPUE likelihood. In the case of the single-stock hypothesis, CPUE1 is overestimated in the ﬁnal years, while CPUE2 is underestimated in the same period. The better ﬁt for the two-stock hypothesis follows through to the improved precision of the estimates of model parameters and management quantities. The trajectories of the exploitable biomass of the H. reedi presents similar patterns in both of the selected models (Figs. 5–7). The single-stock model estimates that between 1989 and 2002 the biomass was less than the BMSY (31,778 t) and from 2001 to 2008 the biomass increased, but recent estimates of biomass are uncertain. Also, this model shows that between 1995 and 1997 the landings level exceeded the MSY (9857 t). Since 1998, the landings have decreased due to catch limits and implementation of closure areas. The two-stock model shows a similar pattern in biomass trajectories. The estimates suggest that the biomass increase started later in the northern zone due to the fact that the reduction in catch occurred later in that zone. The southern zone shows a recent leveling off in biomass, which corresponds to a recent increase in catches in that zone. Furthermore, there is greater uncertainty in the biomass estimations in the southern zone in recent years. This model shows that in the northern zone the landings exceeded the MSY during 1997 and 1999, while in the southern zone, the landings exceeded the MSY during 1995 and 1996.

4. Discussion We ﬁt a total of 44 models for the representation of the biomass dynamics of the shrimp H. reedi (combinations of hypothesis/models/scenarios). From these two structural hypotheses, we select the best model of each one (based on AIC). For both hypothe-

ses the best model ﬁt was the one with six parameters (the single-stock hypothesis that estimates B0 , p, r, , q1 , q2 , and the two-stock hypothesis that estimates B01 , B02 , p, r, , q). For each hypothesis, the best model ﬁt was the one with a symmetric production curve (BMSY /B0 = 0.5). The two-stock model selected was clearly better than the single-stock model selected (single-stock model AIC = −28.57 and two-stock model AIC = −46.58). The lower uncertainty in the biomass estimation in the two-stock model is a consequence of the lower uncertainty in the model’s parameters than the single-stock model (Table 3). In fact, the CVs of the single-stock model’s parameter are almost two times the CVs of the two-stock model’s parameters (with the exception of of CPUE). If the two stocks were assessed independently (i.e. no parameters shared among the two stocks) the CVs were also greater (Table 4).

Table 4 Parameters, management quantities, standard errors, coefﬁcient of variation (CV%) and 95% conﬁdence limits for the two-stock hypothesis with all parameters estimated. Parameter

Value

Std. Err

CV%

Lower bound

Upper bound

B01 B02 p1 p2 r1 r2 1 2 q1 q2 BMSY1 BMSY2 MSY1 MSY2

25,286 33,053 0.3 0.26 0.33 0.33 0.12 0.12 1.53E−05 1.73E−05 12,643 16,526 4213 5469

8422 8007 0.06 0.04 0.09 0.06 0.02 0.02 2.64E−06 2.24E−06 4211 4003 333 468

33% 24% 20% 15% 27% 18% 17% 17% 17% 13% 33% 24% 8% 9%

8779 17,359 0.182 0.182 0.154 0.212 0.081 0.081 1.01E−05 1.29E−05 4389 8680 3560 4552

41,793 48,747 0.418 0.338 0.506 0.448 0.159 0.159 2.05E−05 2.17E−05 20,897 24,372 4866 6386

196

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

Fig. 6. Results from the two-stock hypothesis. (a) Exploitable biomass estimated with the surplus-production model (solid line, with 95% conﬁdence intervals), BMSY (horizontal solid line), northern zone; (b) landings (solid line) and MSY (horizontal dashed line), northern zone. (c) Exploitable biomass estimated with the surplus-production model (solid line, with conﬁdence intervals), BMSY (horizontal solid line), southern zone; (d) landings (solid line) and MSY (horizontal dashed line), southern zone.

These results indicate that having the right spatial structure and sharing information between areas can improve the precision of the estimates of management quantities when multiple indices of abundance are available from different areas and are somewhat inconsistent. The information in the data support values of the shape parameter much greater than the values suggested by the demographic parameters. As mentioned earlier, the demographic parameters support a value of BMSY /B0 closer to 0.3. The AIC was 1.87 units smaller for the selected two-stock model when the shape parameter was estimated (BMSY /B0 = 0.6) than when it was ﬁxed corresponding to BMSY /B0 = 0.3. In addition, the results are substantially different between the two values of BMSY /B0 , and imply different depletion and MSY levels, which would translate into different management actions. The values of the depletion level

Fig. 7. Depletion levels (Bt /BMSY ) for best models of the single- and two-stock hypotheses (Stock Total correspond to the single-stock hypothesis; northern and southern zone correspond to the two-stock hypothesis).

(biomass as a ratio of the virgin biomass, Bt /B0 ) in 1989 are unrealistically low (0.07) for the m = 0.68 model. Often there is not enough information to estimate the shape parameter (Pella and Tomlinson, 1969; Fletcher, 1978; Rivard and Bledsoe, 1978; Hilborn and Walters, 1992; Prager, 2002), which appears to be the case in the H. reedi application, and estimates may be due to quirks in the data rather than information about the shape parameter. Since the best model estimates that the stocks are well above the BMSY levels, increasing catch levels to MSY would reduce the population sizes, providing more contrast in the CPUE time series, which may provide additional information to estimate the model parameters. However, the consequences of different harvest strategies must be carefully analyzed, because catch levels slightly greater than the MSY, in the mid 1990s, probably caused the decline of the population levels. The results of this study identify an inadequacy of the standard hypothesis testing framework. For example, it would be conventional to select the Schaefer model (BMSY /B0 = 0.5) as the null hypothesis because it describes a linear relationship between biomass and the production per unit of biomass. In this case, the analysis of the two-stock model would retain the null hypothesis (i.e. the Schaefer model). However, if a more realistic (based on the arguments of Maunder, 2003) null hypothesis was chosen (BMSY /B0 = 0.3), the null hypothesis would be rejected and the model with BMSY /B0 = 0.6 would be chosen. Therefore, the model used depends on which model was chosen as the null hypothesis. This supports the use of model averaging over hypothesis testing. Standard stock assessment models using a Beverton and Holt (1957) stock–recruitment relationship have a maximum BMSY /B0 = 0.5 (Thompson, 1992), suggesting that the H. reedi stock may have a stock–recruitment relationship that differs from that of Beverton and Holt. The results suggest that the demographic parameters are similar between the two stocks. The only parameter that differed was B0 ,

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

197

Fig. 8. Area of operation (nm2 ) of the ﬂeet ﬁshing for Heterocarpus reedi between 1995 and 2007.

which was greater for the southern zone than for the northern zone. It is possible that B0 would also be the same for the zones if it were put in terms of virgin biomass per unit of suitable habitat. A proxy of the potential habitat availability for the species, the area between 200 and 600 m in each zone, shows larger habitat availability in the southern zone, with an area of 2512 (nm)2 , than the northern zone with an area of 1870 (nm)2 . Based on the results from the two-stock model, this corresponds to 15 and 14 t per nm carrying capacity for the northern and southern zones, respectively, indicating that the virgin biomass parameter could be shared between the two zones if modeled on a per habitat basis. On the other hand, the area of operations of the ﬂeet has had important variations during the history of the ﬁshery. These changes are due mainly to management policies, economic incentives of the industry, and levels of abundance of the species (Fig. 8). It is possible that the two zones may be linked through recruitment. There is no information available about the larval and postlarval stages of H. reedi. However, another species of this genus, Heterocarpus sibogae, has a larval period (before metamorphosis) of between 33 and 41 days (Iwata et al., 1986). Therefore, it is possible that H. reedi has similar larval stage duration. The habitat where H. reedi inhabits has a very dynamic system of currents, with the Humboldt Current’s coastal branch moving toward the north and the coastal Peru-Chile undercurrent moving toward the south (Fuenzalida et al., 2007). It is therefore possible that the two zones may have some degree of mixing due to the advection of larval from north to south and vice versa. Surplus-production models combine all the biological processes into a single function, so it is not possible to separate recruitment and compare the estimates of recruitment between the two zones or to use a

joint stock–recruitment relationship. A more detailed population dynamics model (e.g. an age-structured model) and data that provide information on recruitment (e.g. catch-at-age data) would be needed in investigating the linkage in recruitment between the two zones. It appears that the main difference between the two zones is the exploitation history, with the catch in the southern zone reducing earlier than in the northern zone and consequently the biomass of the southern zone increased earlier than did that of the northern zone. This implies that local depletion can occur in this stock and that differences in management among zones may require implicitly modeling sub-stocks in the assessment of this species. Our results show that sharing information among sub-stocks can increase the precision of the estimates of management quantities. We only used two stocks and it is possible that a ﬁner spatial structure may be more appropriate. A population dynamics model with more ﬁner spatial scale stocks that shares information by treating parameters as random effects may further improve the parameter estimates and management advice. Acknowledgements We thank William Bayliff, Andre Punt and two anonymous reviewers for providing comments that improved the ﬁnal version of the manuscript. Appendix A. Tables of parameter estimates See Tables A1 and A2.

Table A1 Models for one stock hypothesis, biomass dynamic of Heterocarpus reedi. Parameters and statistics

Model 1, m = 2 (BMSY /B0 = 0.5)

Number of parameters/extra parameters 5 Base

6 sd2

6 q2

7 sd2, q2

−ln(L) AIC

−37.81 −27.81

−38.98 −26.98

−40.57 −28.57

−41.49 −27.49

B0 p r 1 2 q1 q2

61,845 0.257 0.320 0.212 0.212 8.4E−06 8.4E−06

37,986 0.267 0.495 0.287 0.141 1.2E−05 1.2E−05

63,556 0.261 0.310 0.197 0.197 7.5E−06 8.8E−06

40,391 0.272 0.464 0.252 0.141 1.0E−05 1.2E−05

198

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

Table A1 (Continued ) Parameters and statistics

Model 2, m = 0.68 (BMSY /B0 = 0.3)

Model 3, m = free parameter

Number of parameters/extra parameters 5 Base

6 sd2

6 q2

7 sd2, q2

BMSY MSY

30,922 9898

18,993 9408

31,778 9857

20,196 9378

−ln(L) AIC

−37.15 −27.15

−37.18 −25.18

−39.87 −27.87

−39.88 −25.88

B0 p r 1 2 q1 q2

5.56E+07 0.000 0.016 0.216 0.216 7.3E−06 7.3E−06

4.18E+25 0.000 2.9E−08 0.203 0.231 6.9E−06 6.9E−06

1.05E+07 0.002 0.028 0.200 0.200 6.5E−06 7.6E−06

460,871 0.038 0.113 0.220 0.181 7.1E−06 8.3E−06

BMSY MSY

16,670,678 265,732

1.3E+25 3.6E+17

3,139,260 89,167

138,092 15,551

−ln(L) AIC

−38.78 −26.78

−39.11 −25.11

−41.53 −27.53

−41.74 −25.74

B0 p r m1 1 2 q1 q2

37,088 0.361 0.429 4.517 0.207 0.207 9.8E−06 9.8E−06

35,494 0.289 0.500 2.535 0.273 0.149 1.2E−05 1.2E−05

39,145 0.363 0.407 4.509 0.191 0.191 8.6E−06 1.0E−05

37,138 0.309 0.463 2.890 0.231 0.153 9.9E−06 1.2E−05

BMSY MSY BMSY /B0

24,157 10,373 0.651

19,363 9691 0.546

25,485 10,376 0.651

21,180 9808 0.570

Bold values means best ﬁts.

Table A2 Models for two-stock hypothesis, biomass dynamic of Heterocarpus reedi. Model 1, m = 2 (BMSY /B0 = 0.5)

Model 2, m = 0.68 (BMSY /B0 = 0.3)

Parameters and Number of parameters/extra parameters statistics 5 6 6 6 r2 Bstart2 Base B02

6 sd2

6 q2

7 B02 , r2

7 7 B02 , Bstart2 B02 , sd2

7 B02 , q2

10 B02 , r2 , p2 , sd2, q2

−ln(L) AIC

−29.98 −19.98

−58.58 −46.58

−53.52 −41.52

−56.44 −44.44

−32.58 −20.58

−47.43 −35.43

−58.74 −44.74

−58.61 −44.61

−58.58 −44.58

−58.88 −44.88

−59.13 −39.13

B01 B02 p1 p2 r1 r2 sd1 sd2 q1 q2

28,349 28,349 0.476 0.476 0.317 0.317 0.264 0.264 1.02E−05 1.02E−05

28,578 35,227 0.254 0.254 0.315 0.315 0.119 0.119 1.63E−05 1.63E−05

35,618 35,618 0.241 0.241 0.263 0.322 0.137 0.137 1.53E−05 1.53E−05

44,480 44,480 0.167 0.230 0.266 0.266 0.126 0.126 1.55E−05 1.55E−05

28,046 28,046 0.482 0.482 0.322 0.322 0.168 0.376 9.16E−06 9.16E−06

23,363 23,363 0.379 0.379 0.395 0.395 0.162 0.162 1.08E−05 1.82E−05

26,545 34,151 0.263 0.263 0.333 0.318 0.119 0.119 1.66E−05 1.66E−05

28,668 34,690 0.250 0.259 0.316 0.316 0.119 0.119 1.64E−05 1.64E−05

28,558 35,205 0.254 0.254 0.315 0.315 0.120 0.118 1.63E−05 1.63E−05

26,919 32,985 0.266 0.266 0.327 0.327 0.118 0.118 1.61E−05 1.70E−05

25,286 33,053 0.298 0.259 0.333 0.331 0.118 0.117 1.53E−05 1.73E−05

BMSY1 BMSY2 MSY1 MSY2

14,174 14,174 4493 4493

14,289 17,613 4500 5547

17,809 17,809 4690 5736

22,240 22,240 5917 5917

14,023 14,023 4509 4509

11,681 11,681 4614 4614

13,272 17,076 4413 5438

14,334 17,345 4535 5488

14,279 17,602 4500 5547

13,460 16,493 4396 5387

12,643 16,526 4213 5469

Parameters and statistics

Number of parameters/extra parameters 5 Base

6 B02

6 r2

6 Bstart2

6 sd2

6 q2

7 B02 , r2

7 B02 , Bstart2

7 B02 , sd2

7 B02 , q2

−ln(L) AIC

−32.04 −22.04

−55.78 −43.78

−51.02 −39.02

−55.21 −43.21

−35.03 −23.03

−47.66 −35.66

−56.08 −42.08

−56.07 −42.07

−55.90 −41.90

−56.07 −42.07

B01 B02 p1 p2 r1 r2 sd1

33,565 33,565 0.358 0.358 0.461 0.461 0.249

1,263,810 1,543,580 0.007 0.007 0.048 0.048 0.129

337,544 337,544 0.031 0.031 0.079 0.095 0.147

62,684,100 62,684,100 0.000 0.000 0.012 0.012 0.131

30,789 30,789 0.361 0.361 0.507 0.507 0.152

38,710 38,710 0.262 0.262 0.378 0.378 0.161

393,439 509,065 0.021 0.021 0.081 0.076 0.128

468,618 532,908 0.018 0.020 0.075 0.075 0.128

1,004,060 1,225,480 0.009 0.009 0.052 0.052 0.121

430,979 521,775 0.020 0.020 0.076 0.076 0.128

C. Montenegro et al. / Fisheries Research 100 (2009) 191–199

199

Table A2 (Continued ) Model 2, m = 0.68 (BMSY /B0 = 0.3)

Model 3, m = free parameter

Parameters and statistics

Number of parameters/extra parameters 5 Base

6 B02

6 r2

6 Bstart2

6 sd2

6 q2

7 B02 , r2

7 B02 , Bstart2

7 B02 , sd2

7 B02 , q2

sd2 q1 q2

0.249 1.10E−05 1.10E−05

0.129 1.40E−05 1.40E−05

0.147 1.30E−05 1.30E−05

0.131 1.43E−05 1.43E−05

0.362 1.05E−05 1.05E−05

0.161 1.01E−05 1.63E−05

0.128 1.44E−05 1.44E−05

0.128 1.43E−05 1.43E−05

0.137 1.39E−05 1.39E−05

0.128 1.37E−05 1.46E−05

BMSY1 BMSY MSY1 MSY2

10,057 10,057 4639 4639

378,680 462,508 18,010 21,997

101,139 101,139 7960 9603

18,782,256 18,782,256 224,288 224,288

9225 9225 4680 4680

11,599 11,599 4387 4387

117,888 152,533 9573 11,659

140,414 159,677 10,508 11,950

300,850 367,195 15,756 19,230

129,136 156,341 9825 11,894

Parameters and Number of parameters/extra parameters statistics 6 7 7 7 r2 Bstart2 Base B02

7 sd2

7 q2

8 B02 , r2

8 8 B02 , Bstart2 B02 , sd2

8 B02 , q2

8 B02 , m2

−ln(L) AIC

−33.26 −21.26

−59.65 −45.65

−54.30 −40.30

−57.27 −43.27

−36.54 −22.54

−48.08 −34.08

−59.68 −43.68

−60.13 −44.13

−59.78 −43.78

−59.77 −43.77

−60.32 −44.32

B01 B02 p1 p2 r1 r2 m1 m2 sd1 sd2 q1 q2

47,800 47,800 0.237 0.237 872,727 872,727 1.5E−07 1.5E−07 0.241 0.241 1.13E−05 1.13E−05

21,501 27,038 0.316 0.316 0.361 0.361 3.246 3.246 0.116 0.116 1.72E−05 1.72E−05

29,873 29,873 0.285 0.285 0.288 0.356 3.347 3.347 0.134 0.134 1.56E−05 1.56E−05

29,788 29,788 0.246 0.324 0.329 0.329 3.881 3.881 0.124 0.124 1.61E−05 1.61E−05

42,116 42,116 0.241 0.241 966 966 1.53E−04 1.53E−04 0.145 0.350 1.12E−05 1.12E−05

28,013 28,013 0.332 0.332 0.391 0.391 1.249 1.249 0.160 0.160 1.06E−05 1.75E−05

20,983 26,814 0.318 0.318 0.368 0.362 3.191 3.191 0.116 0.116 1.74E−05 1.74E−05

18,004 26,373 0.382 0.318 0.381 0.381 4.393 4.393 0.114 0.114 1.72E−05 1.72E−05

21,266 26,790 0.316 0.316 0.365 0.365 3.308 3.308 0.122 0.108 1.74E−05 1.74E−05

21,093 26,417 0.318 0.318 0.366 0.366 3.120 3.120 0.115 0.115 1.71E−05 1.77E−05

21,148 27,056 0.319 0.319 0.369 0.369 2.555 4.259 0.114 0.114 1.70E−05 1.70E−05

BMSY /B01 BMSY /B02 BMSY1 BMSY2 MSY1 MSY2

– – – – – –

0.592 0.592 12,729 16,006 4597 5780

0.598 0.598 17,854 6009 5134 2139

0.625 0.625 18,604 5685 6126 1872

– – – – – –

0.409 0.409 11,470 5997 4488 2346

0.589 0.589 12,355 5452 4544 1973

0.646 0.646 11,639 5545 4431 2111

0.596 0.596 12,664 5477 4620 1998

0.585 0.585 12,332 5409 4509 1978

0.547 0.641 11,569 5577 4274 2060

Bold values means best ﬁts.

References Beverton, R.J.H., Holt, S.J., 1957. On the Dynamics of Exploited Fish Populations. Fish. Invest. Ser. II, vol. 19. Mar. Fish. G.B. Minist. Agric. Fish. Food, 533 pp. Fletcher, R.I., 1978. On the restructuring of the Pella–Tomlinson system. Fish. Bull. 76, 515–521. Fox, W.W., 1970. An exponential surplus-yield model for optimizing exploited ﬁsh populations. Trans. Am. Fish. Soc. 99, 80–88. Fuenzalida, R., Schneider, W., Blanco, J.L., Garcés-Vargas, J., Bravo, L., 2007. Sistema de corrientes Chile-Perú y masas de agua entre Caldera e Isla de Pascua. Ciencia y Tecnología del Mar 30 (2), 5–16. Gilbert, D.J., 1992. A stock production modeling technique for ﬁtting catch histories to stock index data. New Zealand Fisheries Assessment Res. Doc. 92/15. Hilborn, R., Mangel, M., 1997. The Ecological Detective: Confronting Models with Data. Princeton University Press, Princeton, NJ. Hilborn, R., Walters, C.J., 1992. Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty. Chapman and Hall, New York, 570 pp. Iwata, Y., Sugita, H., Kobashi, T., Deguchi, Y., 1986. Larval development of the pandalid shrimp Heterocarpus sibogae De Man. Bull. Coll. Agric. Vet. Med. Nihon Univ. 43, 140–150. McCullagh, P., Nelder, J.A., 1989. Generalized Linear Models. Chapman and Hall, London. Maunder, M.N., 2003. Is it time to discard the Schaefer model from the stock assessment scientiﬁc’s toolbox? Fish. Res. 61, 145–149.

Maunder, M.N., Starr, P.J., 1995. Sensitivity of management reference points to the ratio of BMSY /B0 determined by the Pella–Tomlinson shape parameter ﬁtted to New Zealand rock lobster data. New Zealand Fisheries Assessment Res. Doc. 95/10. Pella, J.J., Tomlinson, P.K., 1969. A generalized stock production model. Inter-Am. Trop. Tuna Comm. Bull. 13, 421–496. Polovina, J.J., 1989. A system of simultaneous dynamic production and forecast models for multispecies or multiarea applications. Can. J. Fish. Aquat. Sci. 46, 961–963. Prager, M.H., 2002. Comparison of logistic and generalized surplus-production models applied to swordﬁsh, Xiphias gladius, in the north Atlantic Ocean. Fish. Res. 58, 41–57. Rivard, D., Bledsoe, L.J., 1978. Parameter estimation for the Pella– Tomlinson stock production model under nonequilibrium conditions. Fish. Bull. 76, 523–534. Roa, R., Ernst, B., 1996. Age structure, annual growth, and variance of size-at-age of the shrimp Heterocarpus reedi. Mar. Ecol. Prog. Ser. 137, 59–70. Schaefer, M.B., 1954. Some aspects of the dynamics of populations important to the management of commercial marine ﬁsheries. IATTC Bull. 1, 25–56. Seber, G.A.F., 1973. The Estimation of Animal Abundance and Related Parameters. Charles Grifﬁn and Company, London, 506 pp. Thompson, G.G., 1992. Management advice from a simple dynamic pool model. Fish. Bull. 90, 552–560.