An individual based model of female brown crab movements in the western English channel: modelling migration behaviour

Abstract

An individual based model (IBM) of the female brown crab Cancer pagurus population exploited off South Devon, UK is described. Size dependent movement rules are ascribed to individuals based on previous observations of predominantly westward migration down the English Channel. Two additional versions of the movement rules explored whether the empirically derived rule was necessary to model the temporal and spatial distribution of crabs. Local crab movement was dependent on substrate type and water depth. Females prefer a soft substrate in which they can bury when temperatures are low or they have eggs to incubate. Crabs have size dependent depth preferences with larger crabs preferring greater depths. Two recruitment functions are used which relate the number of incoming crabs to the sea surface temperature five years earlier. Model outputs were tested against 10 years of logbook data from three crab fishers and against data from a year-long sampling programme on eight of the vessels exploiting the area. The model reproduces the long-term pattern which is mostly temperature driven. Spatial variation in catch is captured effectively by the model with more crabs being caught in the east of the area than the west and more caught offshore than inshore. The significance of the results is discussed in relation to the crab life cycle, management of the fishery and the potential effects of increasing temperatures.

Introduction

Marine animals living in dynamic fluid environments have had to evolve life cycles that allow them to persist in a contained ‘home'area. For example, cod (Gadus morhua, Linnaeus 1758) in the Barents Sea move south against the prevailing current to the Lofoten Islands to spawn. Their eggs and larvae are then transported by the remnants of the North Atlantic Drift back into the Barents Sea (Hylen et al., 2008). Plaice (Pleuronectes platessa, Linnaeus 1758) in the southern North Sea spawn in the Southern Bight and eggs and larvae are carried by the NE flowing current up towards the Wadensee area. Here, the young fish move to shallow water where they grow until they are large enough to move into the central North Sea. Adults then migrate back against the current to the Southern Bight in time to spawn (Cushing, 1990; Hunter et al., 2004). There is some evidence that female brown crabs off the Northeast coast of England move north before they lay their eggs, with the resulting larvae being carried south in the prevailing current (Nichols et al., 1982; Eaton et al., 2003). In addition, there is evidence that both Dungeness crabs (Cancer magister, Dana 1852) (Rasmuson, 2013) and red crabs (Pleuroncodes planipes, Simpson 1850) (Longhurst, 1968; Robinson et al., 2004) have movement behaviours that keep them within a limited area, despite their drifting early larval stages.

A recent study of brown crab (Cancer pagurus, Linnaeus 1758) migration (Hunter et al., 2013) has shown that female crabs move almost exclusively east to west down the English Channel (Figure 1 from Hunter et al.,). This movement is shown for crabs tagged in the eastern Channel and off the Devon coast. Although not explicit from the results of the tagging programme, the migration pattern is similarly suggestive of counter current behaviour, where it can be hypothesised that the females move up current towards the west so that their planktonic larvae can be carried back up the Channel as they develop. When temperatures drop to around 13 °C towards the beginning of December, brooding females bury in the sediment for around 180 days, during which period they remain inactive and do not forage (Hunter et al., 2013).

Figure 1. The Inshore Potting Agreement area (2013) (Trawling is permitted in Area 3 provided vessels have an engine power of no more than 100kw and scallop dredging is permitted provided vessels use no more than two tow bars <2.6 metres in total length and there are no more than three dredges attached to each tow bar.) Zone 1: Trawling 1st Jan – 31st March. Zone 2: Trawling 1st Jan – 1st June. Zone 3 (1): Trawling 1st Jan – 31st March. Zone 3 (2): Trawling 1st Jan – 31st March. Zone 4: Trawling 1st Feb – 31st August. Zone 5: Trawling all year. Corridor Trawling 1st Feb – 31st March. Areas 1, 2 and 3- Pots all year. The grid in red with 8 regions was used for collecting data on board fishing vessels (see text).

Crab movements and the influence of temperature imply that the spatial distribution of crabs in the Channel is likely to be strongly influenced by these two factors. To explore the consequences of migration and temperature for the crab fishery we constructed an individual based model (IBM) of female crabs. This model gives the crabs a set of rules governing their movements and how they should respond to temperature variation through the year. To test the temporal and spatial patterns emerging from the model we used data from the south Devon crab fishery which is carried out from the ports of Salcombe and Dartmouth in an area largely reserved for crab fishers. The work was carried out in cooperation with fishers from the South Devon and Channel Shellfishermen's Association (SDCSA) (Pearson et al., 2020) building on earlier work by Blyth and others (e.g. Hart, 1998; Blyth et al., 2002), to gather data on female crab spatial and temporal distribution. Our work involved a year-long sampling programme on board a selection of fishing vessels along with a collation of data from the logbooks of three fishers to produce a 10-year sequence of catches.

The area exploited by the crabbers is shown in Figure 1. The fishery takes place within the United Kingdom 6-mile territorial limit and an agreement established in 1975 between the potters and fishers using towed gear, defined an area within which only potting can take place. This is known as the Inshore Potting Agreement (IPA). As has been shown by Blyth et al., 2002, 2004; Blyth-Skyrme et al., 2006; Kaiser et al., 2000, 2007, the area closed to towed gear has had significant conservation effects on the habitat and the animals living in it.

The accessibility of fast computers has made possible the development of individual- or agent-based models (IBM/ABM) which have provided new openings for the development of models for studying animal population dynamics (Huse et al., 2002; Grimm and Railsback, 2005; Lindqvist et al., 2020). Most animal populations and the ecosystems of which they are a part, are best approached as complex adaptive systems (CAS) which are characterised by three properties (Holland, 1995; Levin, 1998); sustained diversity and individuality of components, localised interactions among these components and an autonomous process that selects a subset for replication or enhancement from the results of interactions amongst these components. In this paper, we focus on the ability of an IBM to account for the spatial and temporal behaviour of the female crab population in a small area off South Devon, UK and as such the model is not a true CAS. We propose though that the model could form the basis of a fisher directed management system. There is no feedback loop between crab behaviour and the successful achievement of a goal such as reproduction. A future version of the model could incorporate this aspect of a population of individuals so creating a genuine CAS.

Description of the model

Earlier IBMs of crab behaviour have been published by Butler et al. (2001); Butler (2003) for the Caribbean spiny lobster (Panulirus argus, Latreille 1804), by Meynecke and Richards (2014) for the giant mud crab (Scylla serrata, Forskål 1775) and by Bunnell and Miller (2005) for the Chesapeake blue crab (Callinectes sapidus, Rathbun 1896). The first model is implemented in Fortran, the second in R and the third report does not specify the modelling software. In all cases the published papers provide no details of how the model has been constructed making it difficult to replicate and run the model without contacting the authors. Through the work of Grimm and Railsback (2005); Railsback and Grimm (2019) and the ComSES (https://www.comses.net) organisation there is a move to encourage those using IBMs to describe their models using a standard approach and subsequently to make them available on the web so that interested parties can either use or develop the models for their own research.

To this end we describe our model using the Overview, Design concept and Details (ODD protocol) scheme as described by Railsback and Grimm (2019). The model is implemented using the NetLogo programming language (Wilensky, 1999; Wilensky and Rand, 2015). NetLogo has four types of agent; globals, patches, turtles and links. Patches have fixed locations but can have individual properties. Turtles, which can be renamed and in our case are called crabs, can move and also have individual properties. Links are, as the name implies, agents that connect either patches, turtles or both to some central location or to each other. We do not use links in the crab IBM. Global variables are available to all agents.

Overview

Purpose and patterns

The model is designed to replicate a female crab's movement rules, as it progresses through the IPA area off south Devon between the ports of Dartmouth in the east and Salcombe in the west (Figure 1). The aim is to create a model that can be used to test the hypothesis that a crab's movement rules coupled with important environmental variables lead to the temporal and spatial distribution of the crabs as revealed by the catch pattern in the crab fishery. An additional aim is to provide a prototype of a model that could be used to manage the crab fishery.

What patterns are used to determine how useful the model is for its purpose?

We use data collected from the fishery to validate the output from the model. Data gathered on board the fishing boats (Pearson, 2017) were used to examine the relationship between the spatial and temporal pattern of the fishery against the model-generated patterns. Three observations derived from the fishery were compared with the model output:

1. 1. That the highest catch per unit area of small crabs is found in the eastern sector of the IPA (see Figure 1 for IPA area).

2. 2. More crabs were discarded, 0–3 nm offshore compared to 3–6 nm offshore, mostly because they were less than the minimum landing size.

3. 3. The logbook records from three crab fishers were analysed to demonstrate the annual pattern of catches over a ten-year period. This pattern was then compared with the output from the model run over the same time span.

The results of the tagging programme of female crabs (Hunter et al., 2013), were used to devise the principal set of rules for crab movements. Variations of the westerly movement were used to determine whether the observed rules are responsible for the spatial and temporal output from the model.

Entities, state variables and scales

Model Domain

The IPA area and some of the sea area surrounding it is divided into 120 (W-E) × 80 (N-S) patches, each with side 500 m (Figure 2). Land is represented by 4962 patches and sea by 4718 patches. Patches representing sea are labelled by a code characterising whether they are in the IPA. Each patch also has a substrate type and a depth (Supplementary Figures 2 and 3).

Figure 2. The NetLogo interface of the model. The orange shapes represent individual crabs on the seabed. The graphs on the right show the catch and the number of crabs as the model runs through 15 years. The top graph on the right shows the sea surface temperature as it fluctuates whilst the straight line shows the temperature below which the crabs bury and stop moving. On the remaining two right hand graphs the red line gives the total number, blue the number of small crabs, purple medium and green is for large crabs. Tabs on the left provide input on the probability of natural mortality, the probability of the three crab sizes being caught, the temperature at which hibernation begins, the number of years that the simulation should continue for, the number of catch samples that should be taken through a year and the vessels whose catch will be recorded (switch on or off).

Crabs – Only female crabs are modelled as they form 91% of the catch. The crabs are divided into three size groups, small, medium and large, representing elements of the natural population structure. The size structure of the south Devon crab population is given by Sheehy and Prior (2008) where a cubic form of the von Bertalanffy (1938) equation was fitted to carapace width at age for females. The equation that best fitted the growth curve of females was;

$$L_{cw} = L_{\infty \;}( {1-\;e^{{-}k\;\times \;( {T-t_0} ) }\;} ) ^3$$

where L_cw = carapace width; L_∞ = asymptotic carapace width = 203.4 mm for females; k = 0.46, rate at which the curve approaches the asymptote; T = age in years; t₀ = −0.54, age at which length would be zero.

Using this curve, the carapace width was calculated for each age between 3 and 9 years (see Table 3), which includes most of the crabs vulnerable to traps. It is then assumed that the small crabs in the model represent those aged 3, 4, and 5, medium are crabs aged 6 and 7 and finally ages 8 and 9 make up the large size group. These have mean carapace widths of 133.9 mm, 179.9 mm, and 193.8 mm. Using a relationship between carapace width and weight for female crabs around the Isle of Man (Öndes et al., 2017) these widths are translated into weights with small crabs weighing 357 g, medium 852.7, and large 1061.3 g.

Previous crab-tagging experiments (Brown and Bennett, 1980; Bennett and Brown, 1983; Hunter et al., 2013), have shown that female crabs move from east to west down the English Channel. Using this information, the base movement rule for the crabs is designed to replicate this westward migration. The model is also run with two further variations on the movement rules to determine how unique the base rules are for capturing the patterns shown by crab catches.

Fishers are not modelled individually but their activity is represented by their capture gear (hereafter ‘pots’) which are set in fixed locations for each vessel. There are 20 vessels in the model (Clark, 1998) and each patch that is part of a vessel's ‘territory’ contains 16 pots (see Supplementary Figure 1). The number of patches associated with a vessel is derived from the total number of pots each fisher deploys divided by 16. Two of the vessels are excluded from the analysis of the results as they had pots in two different areas and the modelled catches from these could not be separated. It is assumed that pots are emptied daily.

The time scale is one day (24 h) and crabs make just one decision each day. The programme keeps track of days and years. A run of the programme extends for 15 years. During the first five years, recruitment of new crabs for each crab size is constant (Table 1). This assumes that the majority of crabs recruit to the fished stock at age five (Sheehy and Prior, 2008). From year six, recruitment is calculated from one of two temperature dependent recruitment functions (Kirby et al., 2008).

Table 1. The daily input of the three sizes of crab during the first five years of a model's run

This number was varied for each version of the model to make sure that the numbers of crabs in the model remained steady through the 15-year simulation

ND = Recruitment function 2 not run with these two movement rules and so has no scaling factor or start crab numbers.

The first recruitment equation (1) is

$$R = \left( {\displaystyle{{( {SS \times 0.1( {SST-5} ) } ) } \over {365}}} \right) \times sf$$

where R = Recruits entering per day; SS = Spawning stock at day 330; SST-5 = Mean sea surface temperature for days 1–100 off Plymouth five years earlier. This would be when the crabs were larvae in the plankton; sf = an arbitrary scaling factor which matches the number of crabs recruiting with the values used during the first five years of a run.

Dividing by 365 converts the annual recruitment to a daily value.

A second recruitment function is derived from plotting the catch from the fishery, as recorded in the logbooks of 3 fishers, against the sea surface temperature off Plymouth five years earlier when it is most likely that a year class was in the plankton as larvae (Sheehy and Prior, 2008). Doing the same plot and regression for sea surface temperatures four or six years prior to the year of catch yields non-significant relationships. The chosen function does not involve the spawning stock size. The equation (2) is:

$$R = \left({\displaystyle{{( {-16.93 + 24.82\;( {SST-5} ) } ) } \over {365}}} \right)\times sf$$

where R = recruits entering per day; SST-5 = Mean sea surface temperature for days 1–100 off Plymouth five years earlier; sf = an arbitrary scaling factor which matches the number of crabs recruiting with the values used during the first five years of a run.

As a result of the fixed recruitment in the first five years, the number of crabs caught would not be a proper representation of the influence of environmental factors on catch. The catch for the final 10 years of a 15 year run of the model is used for comparison with the ten years of logbook data from three crab boats.

Process overview and scheduling

For the model there are the following entities with the actions they can execute.

Entity action

Crabs move, choose a substrate, respond to the density of other crabs on a patch they might move to, bury in the sediment to incubate eggs, move into the modelled area, die naturally from all causes, be caught in a pot or leave the area.

The order in which these are performed each day is as follows:

Decide whether to stop moving and bury.

Choose the next patch to move to and then move.

Die from natural causes with a probability read from an interface slider (see Figure 2).

New crabs move into the area from the east and south in two versions of the model but in the third their positions are created at random over the sea area.

At day 330, the number of crabs of the three sizes is counted to determine spawning stock size. This point in the year is chosen as it is at the end of the period when the temperature is still above that at which the crabs stop moving. Once they stop no new crabs move into the area but those already there begin to die off so rendering the population less representative of the spawning stock.

If the crab is in a patch with pots then it is caught with a given probability for each crab size read from a slider.

Patches representing the seabed are either in the IPA or not, do or do not have pots, have a depth and a substrate type. The temperature is the same for all patches but varies from day to day. The temperature dataset comes from a buoy moored off Start Bay and the annual pattern is shown in Figure 3. The programme adds a random element to the temperature which is read from a file.

Figure 3. The mean daily temperature over a year plus and minus its standard deviation measured on a buoy in Start Bay, Devon UK (50°17.50'N, 003° 36.97'W). Also shown is one possible value for the temperature at which female crabs stop moving and bury in the substrate.

The model is run with three different movement rules.

Version 1, labelled as MR1 (Movement Rule 1) RF1 (Recruitment Function 1) or MR1 RF2 (Recruitment Function 2). The crabs enter the modelled area along the east or southern border and then choose to move to one of the five patches that is within the 180° western sector of patches surrounding them and which has the depth that the crab size prefers, a soft substrate and the smallest number of other crabs. If there are no neighbouring patches with the right conditions, then the crabs stay put.

Variations on this basic set up were as follows:

Version 2, MR2 RF1. Crabs move into the area along the eastern or southern border of the area modelled and then move to one of the surrounding eight patches with the right depth for the crab size, a soft substrate and the smallest number of other crabs. If there are no neighbouring patches with the right conditions, then the crabs stay put.

Version 3, MR3 RF1. Crabs originate anywhere in the area, placed at random, and they then move at random to one of the surrounding eight patches without regard to any of the environmental variables.

New crabs appearing on the eastern and southern margins of the modelled area simulate reproduction/recruitment.

Fishers empty their pots every day producing a daily record of catch although this periodicity can be changed by recording catch only on a reduced number of days in a year giving weekly or monthly catch numbers.

Design concepts

Railsback and Grimm (2019) list eleven elements under this heading of the ODD protocol. Some are repetitions of what has already been outlined earlier, whilst others do not apply, so have been left out of what follows.

Basic principles

Basic information on brown crab ecology is given by Neal and Wilson (2008). The model uses a stripped-down version of the crab's life history, focusing on the elements for which we have a good knowledge. The immigration rate is a reflection of recruitment of new crabs from the eastern and southern edges of the modelled area. New recruits come from all three sizes of crab. The size of the crab on entry to the modelled area determines how it moves.

• • Small crabs prefer depths less than 30 m.

• • Medium crabs prefer depths greater than 30 m but less than 40 m.

• • Large crabs prefer water over 40 m deep.

These depths are based on the age-based distribution of females as described by Sheehy and Prior (2008). Females stop moving when the sea surface temperature falls below a set value read from a slider on the NetLogo interface.

Emergence

The distribution of the three different sizes of crabs across the IPA and surrounding area and the annual variation in catches are emergent properties of the model. This in turn produces the spatial and temporal pattern of catches across the IPA.

Objectives

What measure of future success of the agents is incorporated in the decision process? The measure of success for a female crab is to survive until the time comes for it to stop moving and bury in the substrate. For the crabs, the choices we model are mostly about choosing the right habitat. Although not modelled explicitly, these choices are those that are assumed to bring about good growth, strong survival and ultimately successful reproduction. We are assuming that habitat choice is positively correlated with fitness (in the Darwinian sense). The model also takes account of the number of crabs on a patch. Favourable patches with the right substrate are likely to attract more crabs than will less favourable patches. The more crabs there are on a patch the more likely it is that the feeding rate will drop, so a mechanism that allows individuals to assess the number of competitors and change their behaviour accordingly was incorporated.

Once the sea temperature falls below a set level, it is assumed that females begin their static period when they bury in the substrate and incubate their eggs (Howard, 1982). In the model this is reflected by the stationary state of the crabs. They are still subject to natural mortality.

Prediction

A crab predicts future fitness through its decision rule to pick the next patch with the best substrate and the preferred depth. The best substrate in the non-reproductive period of the year is assumed to be the one that gives the best feeding opportunities reflected by choosing the candidate patch with the lowest number of other crabs. The high level of competition on the best sites and the need to reduce this competition by moving predicts that the crabs are driven to move from where they are, so generating the continuous movement observed. In two of the model versions (1 and 2) if there are no patches that have the right substrate or depth then the crab does not move as there would be no improvement by doing so.

Sensing

Crabs are able to sense the substrates and depth of the immediate neighbouring patches and how many other crabs will be on the next prospective patch. How they do these tasks is not explicitly modelled, they just ‘know’.

Interaction

There are no explicit interactions between crabs.

Stochasticity

A random number generator is used to determine if a crab will die on a given day. The daily natural mortality rate is set on a slider and is calculated for each crab on each day as ‘mortality (from the slider) ± a number chosen by the programme between 0 and 0.001’. If this number is less than the mortality rate from the slider, the crab dies. The same is true for capture in a pot so that capture probability q_n is ‘catchability from the slider for each crab size ± a number chosen by the programme between 0 and 1’. This introduction of a stochastic element represents the chance nature of death by natural causes and the fact that the crab encountering a pot is to some degree due to chance. The variations around the mean are normally distributed. In addition, the mean sea surface temperature each day has a standard deviation and as the programme runs it adds or subtracts to the mean a number between 0 and the standard deviation.

The number of crabs recruiting to the population is also subject to a random element, as would be expected in nature.

Observation

For each vessel the outputs include day and year, crab catches as a total and by size, total number of crabs in the modelled area and their partition into the three sizes. The data that is output allows the catches to be assigned to specific areas in the IPA as each vessel has a fixed location. Output of the spawning stock size divided into the three crab sizes along with the recruitments derived from this spawning stock size is sent to a separate data file. Each vessel is given an individual switch on the NetLogo interface (Figure 2). Using this, catch from each vessel can be recorded separately. This allows for a sampling regime using different numbers of vessels all recording their catches at the same intervals as selected on the sampling-interval slider. The model interface also shows the values of all the variables input from sliders, the total crab population, the catch and the mean sea surface temperature (see Figure 2).

Details

Initialisation

Each run starts with a predetermined number of crabs randomly distributed through the sea area modelled. The seabed type is a fixed variable for each patch, input from a file, and sea surface temperature has a fixed annual pattern for each run and is also input from a file. The initial number of new crabs of the three sizes entering from the east is fixed but varies with the three versions (see Table 1). This variation, together with the scaling factor in the recruitment function, ensures that the number of crabs remains comparable through all the 15-year simulations. The numbers entering per day from the south are double these numbers to account for the fact that crabs die if they choose to move outside the 120 × 80 frame of the model world. As a result, around half the crabs generated on the southern border can be expected to die as soon as they are introduced.

Input data

Whether a patch is in or out of the IPA or on land is defined by a colour code read from a data file. This file also contains data on seabed type and depth for each patch representing sea and whether there are pots on the patch (File title: ‘Co-ords IPA + pots + vessels + sed2 + depth.txt’). Sea surface temperature recorded daily from a buoy in Start Bay is input from a separate file (File title: ‘Slapton temps daily mean + st dev.txt’, see Figure 3). A further input file contains sea surface temperature off Plymouth for the first 100 days of each year over a ten-year period (File title: ‘SST Plymouth 95-14 M & SD.txt’). This second set of temperature data is used as input to recruitment functions.

Sliders on the NetLogo interface set the values of the natural mortality rate to be used in a run, from which vessel data will be recorded, the number of samples of catch data to be taken each year and the probabilities of capture for each of the three crab sizes. For the simulations described in this paper data are recorded from all the vessels once a month, natural mortality is constant for all runs at 0.00001 per day and the probability of capture for all three crab sizes is set at 0.5 (but see later).

Submodels

At the start of a run the programme activates a sequence of routines to set up patches, crab shapes, crab pots, lists to store the catches from each vessel, lists to store the total population of crabs, the output file for catches, the output file for the spawning stock size and recruitment, the temperature profile in Start Bay and the sea surface temperature off Plymouth for the recruitment function.

A series of sub-routines, called by a ‘to go’ instruction, then implement aspects of the model. The first just sums the area occupied by the pots of each vessel. These data are used to calculate catch per unit area. The remaining sub models are as follows.

Count spawning stock

On day 330 of each year of the run the number of crabs in each size group are summed. The spawning stock size and resulting recruitment is output to a file called ‘Spawning_stock_size.csv’

Burying temperature

This routine looks up the temperature for the day and adds variation to the mean to give the realised temperature for the day.

Crab movement

For all sizes of crab, each individual moves according to which ever movement rule is being implemented (see earlier descriptions of movement rules).

Natural mortality

At the end of each day each crab has a probability of dying from natural causes.

Incubation temperature

This routine looks up the SST temperature for the year, averaged over the first 100 days, and adds some variation to the mean to give the realised temperature for the year.

Generation of new crabs

New crabs appear on the eastern and southern edges of the modelled area in versions MR1 and MR2 but anywhere in version MR3 and have a size distribution reflecting that of a natural population. The number of medium and large crabs recruiting is those that survive from the small and medium cohort of the previous year.

Crab catch

After the crabs have moved the catch sub-model will come into effect. For each crab, if it has landed on a patch containing pots, it is subject to a probability of capture. The sub-model keeps a tally of how many crabs have been caught by all the pots on a patch since the day before and once a maximum has been reached the probability of a crab being caught falls to zero. After the pots have been ‘cleared’, the probability of capture will be reset to its original value and capture can start again. For those vessels that have been selected up to a total of 20, catch data is recorded in an output file (‘Catches_per_vessel_per_day v2.csv’).

Data gathered from the South Devon fishery for model validation

Methods

The monthly catch per unit area fished by pots was extracted from the logbooks of three crab fishers and is shown in Figure 4. This annual pattern is taken as being representative of the fishery as a whole. The spatial distribution of catches was obtained from a year-long sampling programme where trips on eight vessels were made roughly once a month from April to November during the periods when crab catches were significant (Pearson, 2017). The eight vessels were distributed across the IPA shown in Figure 1.

Figure 4a. The mean catch, in kilograms or numbers, for each month over a 10-year period. If the mean weight of caught crabs is 1 kg then the values for each curve are comparable.

Figure 4b. Kilograms/area from fisher logbooks plus and minus the standard error on the mean together with the mean catch/area from the model with MR1 RF1.

Figure 4c. Mean kilograms/area from logbooks plotted with the mean catch/area for the model with MR1 RF1 plus and minus the standard error on the mean.

Results

The data showed that catches in the east were greater than those in the west although only significantly so for small individuals (Table 2). More of these small crabs were discarded in the area 0–3 nm from the coast than in the area 3–6 nm (Table 2). The crabs were mostly discarded because they were below the minimum landing size for females which is 150 mm carapace width. Crabs that had recently moulted and had soft carapaces were most abundant in July, August, and September.

Table 2. Results of two-way analysis of variance of the total catch per unit area either in weight or numbers and the catch of small crabs per unit area>

The Effect Summary, expressed as -log₁₀(p), provides a clearer picture of how significant the effect of the variable is when probabilities of significance are very small. Both variables were analysed together so the Effect Summary applies to the combined analysis.

Running the model

Determining the number of replicates required and statistical analysis

Because of the variance introduced into the model through a random addition to many of the variables, it was necessary to replicate each treatment. To determine the acceptable minimum number of replications, the model using MR1 RF1 was run 20 times and the standard error on the mean catch calculated for each number of runs from 1 to 20. The standard error on mean catch was reduced by 71% after 10 replications and by 79% after twenty replications. As a result of this analysis it was decided that 10 replicates per condition would be acceptable.

All results of both catch data from fishing vessels and simulations from the model were analysed using the JMP Statistical Package (Sall et al., 2005).

The simulated annual pattern of catches

Ten replicates of the model were run for each of the three movement rules. These were compared to the 10 years of logbook catch data. The latter is in the form of mean catch in weight per month and the weights have been divided by the number of patches in the model that contain pots belonging to each of the three vessels. This produces a standardised value comparable with the total number of crabs caught per unit area fished derived from the model. Model output is in terms of numbers. As described earlier the estimated weights of the three crab size groups has an average of 0.76 kg but it was noted that the Devon crabs are bigger than those off the Isle of Man so it is not unreasonable to assume that the average weight of crabs caught off Devon was close to one kilogram, As a result the number caught is assumed to be equivalent to kilogrammes caught. The logbook data are compared with the model output for each of the three movement rules. In addition, under movement rule 1, a comparison was made between the logbook data and model output using each of the two recruitment functions.

The basic annual pattern of increasing catches as spring advances with a decline in late autumn and into the winter period is shown in all the simulated runs (Figure 4). However, the mid-summer dip in catches shown in the logbook data was absent from the simulations, and no simulations demonstrated the same pattern of variation in peak catches shown by the logbook data through the 10-year period. The mean annual catch per area for each version of the model, plotted against the sea surface temperature over the first 100 days of the year off Plymouth five years earlier, is only significant for the version of the model run with MR1 and RF2 (r = 0.74, df = 9, ANOVA F = 9.84; df = 1 9; p = 0.014).

A plot of catch per area for all four simulations (Figure 4) shows that the runs with MR1 and MR2 have values that match the logbook data when it is assumed that the average weight caught in the simulated runs is 1 kg, whilst the catch per area from the simulation with MR3 is much lower. However, catch per area for all four simulations is correlated with their equivalent from the logbooks and explains a significant amount of the variance at p < 0.0001 (MR1 RF1 r ² = 0.71, MR1 RF2 r ² = 0.70, MR2 RF1 r ² = 0.64, MR3 RF1 r ² = 0.67, df = 119 and p < 0.0001 in all cases).

A Principal Component Analysis of the total catch per area per month in the four simulations and the logbook data shows that the first principal component accounts for 91.2% of the variance. This principal component is strongly correlated with sea surface temperature as measured off Slapton with r = 0.89, N = 119. The regression equation is PC1 = −9.67 + 0.76 Mean temperature. RF2 is based on the relationship between crab catch and the sea surface temperature off Plymouth five years earlier. As a result RF2 is not independent of temperature. Consequently, the catch data from MR1 RF2 was removed from the data table and the Principal Component analysis redone. This resulted in the same degree of correlation as for the previous analysis with r = 0.89, N = 119. The regression equation is PC1 = −8.64 + 0.68 Mean temperature. An ANOVA gave F = 457.5, df = 1, 119 and p < 0.0001 and PC1 accounts for 93% of the variance.

The second principal component is correlated with the logbook data mean catch of the three vessels with Mean (Kg/area total) Logbooks = 3.94 + 2.59PC2 (F = 24.85, df = 1, 119, p < 0.0001). The trend for the simulated catches is the opposite of this so that high catches are correlated with positive values of PC2 but in all three cases, the correlation is not significant.

Spatial analysis of catches

In the model there are 7 vessels in the western area and 11 in the eastern part, defined as being west or east of a line going south from Start Point (see Figure 1). Similarly, there are 11 vessels in the inshore region (0–3 nm) and 7 offshore (3–6 nm). An analysis of the mean number of patches per vessel in the west and east and inshore and offshore shows that there is an average of 93 patches per vessel with pots in the west and only 55 in the east. Inshore the average coverage is 57 patches per vessel and offshore 90.

The spatial catches were analysed in this context of vessel and pot distributions. The statistical analysis of the four versions of the model defined by movement rule and recruitment function is shown in Table 2. For the MR1, RF1 there were fewer small crabs caught per unit area in the west than the east and fewer caught inshore than offshore. The same pattern held for total catch per unit area. For MR1 RF2 the same pattern was displayed between catches per area in the west vs the east and inshore vs offshore. This pattern was not changed for catches per area when MR2 RF1 applied. For the third version of the movement rules, MR3 RF1, there was no difference between west to east catches per area of either small crabs or all crabs. Additionally, inshore catches were less than offshore for both small crabs per area and all crabs per area.

Discussion

The correspondence between fisher recorded and simulated catches per month (Figure 4) and the correlation of the first Principal Component with sea surface temperature supports the conclusion that temperature variation is largely responsible for the annual pattern of crab catches in the South Devon fishery. The decreasing temperature in the early winter period prompts the female crabs to bury themselves in the substrate and to stop moving. As such they become unavailable to the pots.

The distribution of catch along the West–East axis is clear in the onboard catch samples (Table 2) and this distribution is echoed by the simulated results for MR1 and MR2. Version MR3 does not have an inbuilt east to west behaviour and does not show the pattern observed in the data derived from fishing boats. Both MR1 and MR2 have crabs entering only from the east and south which leads to a higher density of crabs in the east. This might also explain why the fishing fleet has more vessels in the east than the west (11 vs 7). Fishers report that when there were fewer vessels exploiting the stock, pots were moved around to follow the changing abundance of crabs across the area. The present density of vessels means that those with established areas in the west are unable to move their pots and must cope with fewer crabs per pot.

In the fishery the number of crabs that are discarded because they are below the minimum landing size of 150 mm carapace width is greater in the inshore than the offshore areas (Table 2). This pattern is not echoed by the simulated catches (Table 2). A possible reason for this difference is that the model does not capture the way in which very small crabs move. As crabs metamorphose into the adult form they drop to the bottom and move closer towards the shore where they can be found under rocks and in crevices right up to the intertidal zone (Sheehy and Prior, 2008). Growing crabs begin to move out into deeper water and the model does not capture this phase of the female crab's life history. It would not be difficult to create a subroutine in the model that dealt with juveniles although there is no data on the age-related sequence of behaviours that these crabs would deploy as they moved from shallow to deeper water.

At present the crab fishery is not subject to a catch quota and the only management measure is the 150 mm minimum landing size for female crabs. To render the model fully useful for management purposes, it would be necessary to represent each age group separately and to have an estimate of the probability of capture (q_n) for each age class. This is not a change that would take a great deal of extra programming effort but the run time for a simulation would increase significantly. With all age groups represented it would be easier to estimate the effects of changing the minimum landing size and evaluating quotas if that became necessary. A limit to this approach would be the present lack of a quick method for ageing crabs.

Even in its present form, the value of q_n for the smallest size group can be changed to represent an increase in the minimum carapace width. The carapace width at age is shown in Table 3 as calculated from the von Bertalanffy equation given by Sheehy and Prior (2008) and reproduced earlier. Also shown in the table, is the positive or negative deviation from each age's carapace width from the minimum landing size. A second column shows this deviation for a minimum landing size raised to 160 mm. At this higher value the deviations at age for ages 3, 4, and 5, which are the age classes that we have assigned to the small crab category, increase by a factor of 1.6. Alternatively, the deviation from the present minimum carapace width for the three age groups is 0.6 of the values derived from the increased carapace width. Throughout the model runs presented in this paper, the probability of capture in a trap for all three size classes of crab has been set at q_n = 0.5. Now 0.6 × 0.5 = 0.3, so changing the q_n value to 0.3 can be used as a proxy for altering the minimum landing size from 150 mm to 160 mm.

Table 3. Carapace width at age calculated from the von Bertalanffy equation given in the text

The table also shows the difference between the carapace width at age from the minimum landing carapace width permissible as it is set now and assuming it was raised to 160mm.

The change in catch per unit area for small crabs when q₁ is changed from 0.5 to 0.3 is shown in Supplementary Figure 4. This demonstrates that an increase in the minimum carapace width allowed would reduce catches of small crabs and the percentage reduction could be estimated from the model.

The model can also be used to design a sampling programme to estimate crab stock size. As described, data can be recorded from a subset of the crabbing fleet by switching the vessel on or off on the NetLogo interface (Figure 2). Using this and knowing the number of crabs in the system it would be possible to choose the minimum set of vessels that would yield a good estimate of crab population density in the area.

The principal weakness in the model at present, in the context of management, is the lack of a well-established recruitment function. We have generated two versions, both of which are predicated on a relationship between recruitment and temperature during the larval phase. However rather than establishing the recruitment functions based on a validated relationship between larval production and the recruitment of animals into the fished part of the stock, these functions had to be adjusted arbitrarily to generate a sequence of long-term catches which matches that observed.

The inclusion of the whole female life cycle would also further extend the model as was done by Bunnell and Miller (2005) for the Chesapeake blue crab fishery. To extend the model would require more information on the timings of life history events such as the movement of late larvae out of the plankton and into the near shore area. As already mentioned, we would also need information on how the young crabs move from the nearshore area to deeper water. An additional element could be the inclusion of males although once again, the knowledge of their behaviour is not well established.

Without feedback between the crab's behaviour and the success or failure at a prescribed task, the current model is not a true CAS. A model that included all aspects of the life history would link movement decisions to reproductive output and over time this could change the movement behaviour to increase fitness. Railsback and Harvey (2020) have provided a method called ‘state and prediction-based theory or SPT’ which allows one to build an evaluation of fitness into an IBM, so exploring the consequences of a set of behavioural options. The method is based on the approach used in stochastic dynamic programming (Clark and Mangel, 2000) which defines an optimal outcome for a behaviour and then determines the sequence of decisions needed to get to the desired end point. In an IBM it is not possible to know what the optimum outcome of a sequence of decisions will be, so Railsback and Harvey (2020) compare the outcomes of a set of feasible decision sequences through a suitable time span and choose the sequence that is predicted to give the greatest final fitness value. This process is repeated at the beginning of each simulated time period so the best decision is continuously updated as both the state of the agent and the environment changes.

The lack of detailed knowledge of ontogenetic behaviour patterns of crabs could be overcome by testing variations of the behaviour patterns in the model and matching the output with observed data. Through this it would be possible to establish a set of hypotheses about the behaviours that could then be used to direct field research.

Finally, the most obvious use of the model as it stands is to explore the implications for the crabs of warming seas. The strong correlation between the variation in annual catch and sea surface temperature, and the response of incubating crabs to changes in winter temperatures, indicate that sea temperature is a major determinant of crab behaviour. It would be easy to change temperature in the model and observe the changes that would take place.