Draft A Lagrangian approach to model movement of migratory species

3 We introduce a Lagrangian movement model that can be used to characterize cyclic migrations of iteroparous fish 4 populations. We demonstrate how movement parameters can be estimated using conventionally available catch at 5 age data and provide a description of the potential bias that may arise from model misspecification, data aggregation 6 and non standardized sampling effort. The model can be extended to incorporate covariates representing biological 7 and environmental forces that alter the distribution and migration range of exploited populations. We expect that this 8 movement model will be a useful tool to model fish migration, to illustrate how fisheries dynamics are affected by 9 fish migration and to be used as an operating model in closed loop simulations. We expect that this movement model 10 will be a useful tool to model fish migration. The model can be used to illustrate how fisheries dynamics are affected 11 by fish migration and could also be used as an operating model in closed loop simulations to test the robustness of 12 management frameworks to spatial structure and connectivity. 13

Two of the main drivers of cyclic migration in marine fish species are seasonal availability of food and spawning behavior (Walters and Martell 2004).In many iteroparous species, these drivers are responsible for a continuous migration cycle between feeding and spawning areas (e.g., Hunter et al. 2003;Costa et al. 2012;Merten et al. 2016).Migration modeling is an important component of fisheries science due to the common presence of migratory behavior in exploited species.In this study, we present a novel method for modeling the migration of fisheries resources and the associated fisheries dynamics that can arise from fish movement.
The migration cycle observed in many iteroparous species can vary in extent (i.e.distance or timing) for subgroups within a population, e.g. for different age or size groups (Ressler et al. 2007(Ressler et al. (2008)), sex (Okamura et al. 2014) or sub stocks (Carlson et al. 2014).This variability can lead to spatial segregation of subgroups within a fish population.When spatial segregation is present, subgroups are susceptible to distinct environmental and ecological drivers leading to differences in natural mortality, recruitment and to additional variability in spatial distribution (Ciannelli et al. 2008).In addition, spatial segregation within a population can cause spatial differences in vulnerability and fishing effort or mortality that are independent from fishing gear or fishing techniques (Martell and Stewart 2014).If ignored, migratory movement and spatial segregation can lead to strong bias in surveys and stock assessment which could result in unreliable management advice (McAllister 1998;Waterhouse et al. 2014).It can also lead to failure of management strategies, particularly when considering space/time closures (Martell et al. 2000;Grüss et al. 2011).The evaluation of impacts of fish movement on fisheries management often requires the use of migration models and simulation studies.See Kerr and Goethel (2014) for a comprehensive review on the application of simulation studies to evaluate the impacts of movement and migration on fish population dynamics and fisheries management.
Migration models are diverse in the fisheries literature (Goethel et al. 2011) and can be grouped based on two numerical methods used to implement them: Eulerian and Lagrangian models.These terms are originally applied to fluid dynamics but have been used to describe migration modeling D r a f t for aquatic resources (Kerr and Goethel 2014;Walters and Martell 2004).The two approaches differ in the way that movement is measured.In the Eulerian approach space is divided in predetermined areas and the movement rate of individuals across such areas is the variable of interest.
Eulerian models have been coupled with stock assessment models and tagging studies (Carruthers et al. 2011;Methot and Wetzel 2013;Sippel et al. 2015) and are useful when there is interest in tracking the net flux of individual across predetermined spatial areas.In the Lagrangian approach, the movement of individuals is tracked through time and space and the movement tracks are the variable of interest (Walters et al. 1999).Lagrangian models rely on the assumption that an underlying force drives movement.For example, this force could be environmental drivers forcing ichthyoplankton dispersal (Lett et al. 2008) or homing behavior driving salmonid runs (Cave and Gazey 1994;Branch and Hilborn 2010).Lagrangian models are not commonly applied to iteroparous species, despite the suitability of the approach to model migratory behavior.Migration hypotheses, such as the cyclic movement between feeding and spawning areas, can be used as the underlying pattern that drives the movement of individuals through space in a Lagrangian model.
In this paper, we describe a Lagrangian movement model designed for iteroparous species that perform cyclic migrations throughout their lifetime.The objective of this model is to provide a way of formalizing movement hypotheses into mathematical models.The resulting model can be used to summarize data and test the validity of alternative migration hypotheses and to represent complex movement dynamics as an operating model in closed loop simulation studies.The model is applied to the Pacific hake (Merluccius productus) as a study case.In particular, we focus on the offshore Pacific hake stock that inhabits the California current system.This stock is believed to perform annual migration cycles between the spawning area off southern California and feeding grounds along the West coast of North America, from Oregon to Southeast Alaska (Ressler et al. 2007(Ressler et al. (2008)).The migration cycle of Pacific hake is partially influenced by the age/length of individuals, with older/larger individuals reaching waters further north (Methot and Dorn 1995;Ressler et al. 2007Ressler et al. (2008)).The migration cycle is a key component to consider for management strategies for this stock.The agreement between Canada and the US defines an aggregate (i.e.non-spatial) harvest control rule, but, given that prosecution of the fishery itself affects the mean D r a f t age/size of fish, it also affects the distribution of the stock and hence, the distribution of the fishery's benefits between the parties.In addition to the model description, we provide a description of data requirements to estimate the model movement parameters.We also show how the model can be extended to incorporate covariates representing biological and environmental factors that alter the distribution and migration range of the populations being modeled.

Methods 3 Movement model framework
We decompose the model into the following sections: population dynamics, movement, and fisheries dynamics.This division is somewhat arbitrary as all three parts are interdependent, but the division is done to ease the description of the model.The time dimension is denoted at t.We assumed monthly time steps, so all quantities within the model were computed twelve times within a year/migrations cycle, but any step length could be used.The indexing of time is cyclic, so that the range {t 0 , ...,t max } is repeated every year.The first month of the migration cycle is also indexed with the year counters (y), so that y − 1 = t 0 in the previous year.The space dimension is considered in three scales: area and fishing ground and territory.The variable area (r) refer to small and equally sized intervals of space, which are used to discretize the variables of interest (biomass, fishing effort, etc.).The area range denotes the limits of the modeled space.Fishing grounds and territories are larger than areas, but are contained within the area range (i.e.within the limits of the modeled space).The fishing grounds correspond to areas where a fleet is known to operate and territories correspond territorial waters of a nation or state.Each territory may contain one or more fishing grounds and the biomass within one territory D r a f t is only accessible to the fleets operating within that territory.Lastly, The group dimension is used to represent parcels within an age cohort that move at different speeds.
In the population dynamics section (Table 3) the processes relating to recruitment, survival, aging and growth are described.These processes are modeled for each group in the population through time and space.The spatial section (Table 4) encompasses the description of movement of individuals at age and group through time .And finally, the fishing dynamics section (Table 5) describes the model effort distribution as a function of spatio-temporal effort scalers and of the spatial distribution of fish biomass (summed over all ages and groups).The spatial distribution of fishing effort is used to generate of spatially explicit fishing mortality.In the movement section, groups are modeled in two ways yielding two versions of the model: single group or multiple groups versions.All model variables are defined on Table 2.
The population dynamics is composed of survival and recruitment processes (Table 3).The model assumes age-specific vulnerability and fecundity as well as Beverton-Holt recruitment occurring at age 1. Recruitment and aging are assumed to happen at the first time step within a migration cycle (Equation T3.2 -cases i-iii) and survival to natural and fishing mortality are calculated at monthly time steps (Equation T3.2 -cases ii-iv).Because of the continuous nature of the model, sometimes some portion of the population might be located outside the boundaries of modeled area and therefore outside the fishing grounds.When this happens, we assumed that individuals are subject to natural mortality only (Equation T3.2 -second term on cases ii-iv).Spawning biomass is combined over all areas and groups (Equation T3.3).
The movement section (Table 4) assumes the individuals in a population are distributed along a unidimensional gradient X and that they perform annual cyclic migrations between spawning and feeding areas.This cyclic migration is modeled with a sine function in which the position of individuals change as a function of time (Figure 1).We have developed two alternative versions: single group and multiple groups.In both versions, the cyclic movement of individuals between spawning and feeding areas is modeled with a sine curve (Equation T4.1).The magnitude of the movement is determined by two parameters representing maximum and minimum positions of the migration cycle (X max,a,g and X min ) and the starting time step of the cycle is given by the D r a f t parameter t 0 .The maximum and minimum positions of the cyclic movement, X max and X min , can be modeled as a function of covariates, such as age, size or environmental drivers.Here, model the maximum position as a logistic function of age (Equation T4.2) and to keep the minimum position was constant for all ages.Once the maximum and minimum positions are defined, the sine curve is used to calculate the mean position of individuals in each group of age a, in time step t ( Xa,t ).The proportion of individuals in each area, i.e small intervals of space, is given by the cumulative normal distribution function given the mean position and standard deviation of the cohort (Equation T4.5) or group (Equation T4.5).
The differences between single and multiple groups versions are in how the individuals are distributed around the mean position ( Xa,t ).When a single group is present, the individuals are assumed to be normally distributed around the mean with variance σ 2 X a (Equation T4.4).In the multiple groups version, each group's position follows a group specific normal distribution and the mean and variance of each group is a function of the overall mean and variance (Equations T4.6-T4.9).The main difference between the two versions is that in the single group version, the distribution of fish regenerates to a normal distribution at every time step, but with a new average position.Although individual fish might experience different fishing mortality depending on location, the regeneration assumption does not allow localized depletion to occur.Instead, fish are redistributed across their range in each time step following a normal density function.When multiple groups are considered, each group is distributed according to a group specific distribution that is much narrower than the overall distribution of individuals at an age class.This mechanism allows certain groups to experience different fishing pressures depending on where the group is located, which might lead to higher or lower fishing pressure over extended periods of time.Because the distribution range of each group is narrower the regeneration only occurs within a narrow range, which allows local depletion to become apparent.
In the fishing dynamics section (Table 5), spatial fishing effort allocation is done with a gravity model (Caddy 1975) .These models assume that effort in each area is a function of the latent yearly and monthly effort in a fishing territory and the attractiveness index of that area (Equation T5.3).
The effort potential quantities E y,k and E m,k are scaling matrices of dimensions Y x K and t max x K D r a f t respectively.The values in these matrices range between 0 and 1.Values of 0 indicate that no effort is allowed in that particular year or month, conversely, values of 1 indicate that full effort is allowed in that particular year or month.The attractiveness index can include factors such as fishing cost, target fish abundance and bycatch abundance (Caddy 1975;Walters and Martell 2004).We make attractiveness (Equation T5.2) a function of vulnerable biomass V B r,t−1 , the power parameter λ and effort potential E pot,r (Equation T5.1).We use the vulnerable biomass in the previous time step because we assume that effort distribution is guided by the biomass distribution observed in the previous time step.The parameter λ is used to indicate if the attractiveness is directly proportional to abundance (λ = 1), or if the fleet tends to disproportionately aggregate in high abundance areas (λ > 1).One example in which λ > 1, is when there is communication between fishing vessels when a high abundance area is located.In such cases, fishing effort would tend to aggregate in high abundance areas generating a patchy distribution of effort.The effort potential E pot,r parameter represents the avoidance factors for a given area, e.g fishing costs and bycatch avoidance.This parameter can be used to represent a range of differences in fishing fleets, such as storage capacity, autonomy at sea and distance between fishing grounds and home port.The E pot,r parameter can also be used to represent areas that are avoided due to high bycatch occurrence.Avoidance areas affect the ability to concentrate fishing effort at a given location and therefore should be considered in the modeling process.Effort is then multiplied by q, the effort scaler, resulting in the fishing mortality in that area.Lastly, catches are calculated for each time step using the Baranov catch equation (Equation T5.4).
Process and observation random error are incorporated in the model.The process random error was represented annual recruitment variability, annual variability in the maximum average position, and annual variability in the effort scaler q.All these variability components were modeled with lognormally distributed error.The age composition data, in numbers and aggregated by fishing ground, are generated with multivariate logistic sampling error.

Application to Pacific hake and simulation-estimation procedure
The offshore Pacific hake stock was used as an example to illustrate the model dynamics.In the Pacific hake case, the fish are known to perform annual migration between the waters off South California and northern British Columbia.Therefore we model the movement of hake in terms of Latitude degrees, from 30 • N to 60 • N (Figure 2).The driving population dynamics parameters for Pacific hake were obtained from the 2015 stock assessment (Taylor et al. 2015).The parameters for the effort dynamics were set to mimic the fisheries dynamics described in the stock assessment document (Taylor et al. 2015).The parameters used in the simulation-estimation procedure are listed in Table 6.
We did a simulation experiment to demonstrate the estimability of the movement parameters.
We simulated and estimated population movement dynamics using both single group and multiple groups versions (20 groups).We simulated total catch and catch at age data with observation error and used that data to estimate the models parameters.
The model parameterization is divided into two categories: fixed (extracted from other models), and parameters that could be estimated, given seasonal catch at age data.The fixed parameters include all the population dynamics and fisheries capacity parameters.These parameters were the recruitment function parameters (R o and h) and natural mortality (M).These parameters were considered as a fixed input to the model and were assumed to be known without error in both simulation and estimation models.The estimable parameters are: t 0 , CV , a 50 , σ X max and q.These parameters were estimated with a multivariate logistic likelihood function fitted to simulated age composition data.
A total of 12 simulation-evaluation scenarios were considered in this study (Table 7).We considered two data aggregation scenarios with data reported from three or five large fishing grounds (aggregated over all areas within fishing ground).These two levels of data aggregation were considered in order to explore the sensitivity of the model to levels of spatial aggregation of the age composition information.We initially considered these two levels of data aggregation for the single and multiple groups version, however we got very low levels of model convergence (50% or less) when 3 fishing ground were considered using the multiple groups version.This fact lead us D r a f t to drop the 3 fishing ground scenarios when the multiple groups version was used.The aggregated catch at age data was assumed to have multivariate logistic error with two levels of observation error, i.e. the standard deviation (τ) being either (0.4 or 1.0).These two levels of variability in the age composition data were chosen in order to evaluate the sensitivity of the model to measurement and sampling error.In addition, we considered two attractiveness scenarios (λ = 1 and λ = 2).
The different levels of attractiveness are likely to change the spatial distribution of fishing effort and would likely impact the degree of depletion experienced by fish in different areas.This is particularly relevant for the multiple group scenarios as the different λ values might affect the local depletion patterns.
The estimation model was identical to the simulation model and parameters were estimated with a multivariate logistic likelihood function fitted to simulated age composition data and a lognormal likelihood fitted to the simulated total catches.A total of 100 simulations were run for each scenario-version combination.Simulation and estimation routines were performed using ADMB (Fournier et al. 2012).The code and simulated data are available for download from https://github.com/catarinawor/Lagrangian 4 Results

Model Dynamics
Figure 3 shows spatial distribution of biomass in the absence of fishing for the single and multiple groups versions.If fishing is absent, the biomass spatial distribution is practically identical.However, the spatial distribution for the two models tends to change if fishing is present and the effort is not homogeneously distributed throughout the natural distribution of the resource (Figure 4).For the multiple groups version, the spatial distribution of fish at each age tends to deviate from the initial normal distribution assumption as the fish grow older.This distortion is caused by the fact that not all groups are subject to the same effort intensity, hence they encounter different fishing mortality, and therefore depletion levels over their lifetime.

D r a f t
In the Pacific hake example, effort is assumed to be concentrated towards the northern areas (bars on Figure 4).Therefore, when the multiple groups version of the model is considered (dark line on Figure 4), the groups that move to higher latitudes tend to be subject to stronger fishing pressure and therefore become more depleted than groups that remain in lower latitudes.This effect cannot be detected for the younger ages (Figure 4, age 1).However, as the cohort ages, the groups that move further (located to the right of the mean distribution for the entire cohort) start exhibiting perceptibly higher depletion levels compared to groups within the same cohort that do not migrate as far (located to the right of the mean distribution for the entire cohort, Figure 4, age 5).The higher depletion levels on the fish that move further also causes the mean in the overall distribution of each cohort to shift to the south, which, overtime would tend to diminish the availability of fish available to the fishing fleets in the northern areas.

Simulation-estimation
The simulation-evaluation analysis showed that the five key parameters of the movement model can be estimated given spatial catch at age data.When data were simulated using the single group version (scenarios 1-8 -Table 7), it was possible to predict the main movement parameters with practically no bias, i.e. parameter estimates were within 10% of true values (Figure 5).Variability in parameter estimates was lower when data were reported for five fishing grounds when compared to scenarios where data were reported only for three fishing grounds.When the data reporting was assumed to occur for five fishing grounds the data was less aggregated and therefore more informative to the estimation of the movement pattern along the latitudinal migration route of the fish.As expected, higher values of τ, the standard deviation for the catch at age error, resulted in higher variability in the parameter estimates.No difference was observed for scenarios with different values for the λ parameter.The parameter λ reflects the capability effort aggregation in areas of high abundance.It is likely that the effects different values of λ were not noticeable because of the regeneration assumption that is inherent to the single group version of the model.In other words, even though effort might have aggregated in certain areas, it did not affect the overall distribution of fish for each cohort.

D r a f t
When data were simulated with the multiple groups version (scenarios 9-12 -Table 7), bias in the estimated parameters became more prominent (Figure 5), although they are relatively low, up to 30% median relative error.The effort scaler parameter, q, was underestimated in all scenarios.
The parameters related to the maximum position at age, a 50 and σ X max were also underestimated.
In order to evaluate the impact of the parameter bias on the predicted biomass distribution, we used plotted the median predicted biomass for scenario 15, one of the cases where the bias was most prominent (Figure 6).Despite the bias in parameter estimates for the multiple groups version the impact on the distribution of total biomass over time predicted by the model was relatively small (Figure 6).Because of the lower effort predictions due to the underestimation of q, higher median biomass distributions are predicted by estimation model.However, little impact is seen in the overall proportion of biomass in each area.

Discussion
The Lagrangian movement model described in this paper is an alternative to traditional Eulerian approaches commonly used to model the distribution of adult iteroparous fish (Goethel et al. 2011;Sippel et al. 2015).The Lagrangian approach shown here allows for the explicit consideration of migration hypotheses, such as the cyclic migration between spawning and feeding grounds, and exploration of potential impacts of covariates in shaping migration variability within a population.
Traditional Eulerian models are generally represented by spatially discrete box or bulk transfer models (Carruthers et al. 2015;Methot and Dorn 1995)  here is continuous in space and therefore does not require the definition of homogeneous spatial areas.However, the model output can be aggregated in larger areas for comparison with historical data or for investigation of management questions relating to political boundaries.
Advection-diffusion models are continuous in space but are much more data intensive and generally require the availability of tagging and tracking data (Sibert et al. 1999;Costa et al. 2012).
Similarly to advection-diffusion models, the Lagrangian model presented in this paper is continuous in space and time, allowing for predictions of biomass at any location in the time-space continuum.However, differently from the diffusion models, the approach shown here does not assume that animals move at random.Instead, the movement is assumed to be directed by an innate migration hypothesis, frequently derived from observed seasonal size or age distribution of the species of interest (e.g., Ressler et al. 2007Ressler et al. (2008;;Francis and Clark 1998).The use of a migration assumption replaces the need to directly estimate advective terms to explain fish movement and, therefore, does not rely on tagging data to determine the direction of fish movement.This is an advantage because for many fish species tagging data is not readily available and tagging studies can be difficult to carry out.This is particularly true for deep water species due the complications and high mortality rates resulting from the barotrauma caused by bringing the individuals to the surface for tagging (Nichol and Chilton 2006;Winter et al. 2007).In addition, tagging studies usually involve high operational costs which include tag deployment, recapture surveys and/or publicity campaigns to increase voluntary tag reporting rates in the fisheries and, in some cases, high costs associated with the tags themselves (Pine et al. 2003) however, that the approach described here is not meant to be a complete replacement for tagging studies.There several caveats associated with the present approach, including the assumption of known population and fishing dynamic parameters.Such parameters are seldom known with certainty and if the assumed values deviate significantly from the reality, the model outputs would likely be severely biased as well.
Appreciable differences between the two versions of the model, i.e single group and multiple groups versions, were observed.These differences are associated with the way in which movement of each age class is treated.In the single group version, all individuals within a cohort are considered equal and individuals are assumed to be normally distributed around the mean position at time.In the multiple groups version each cohort is sliced into groups which have a group-specific mean position at time.The single group model is simpler and, for this reason, much more computationally efficient.However, it relies on the assumption of regeneration of spatial distribution at each time step.That is, fish of a given age are assumed to spread spatially following a normal distribution at each time step regardless of possible local depletion due to concentration of fishing effort.This might be a reasonable assumption if the distribution of effort is relatively homogeneous over the distribution of the age group being exploited.However, effort concentration and local depletion is not uncommon and has been observed for many fish stocks around the world (e.g., Maury and Gascuel 2001;Bez et al. 2006).Alternatively, if fishing effort is known to be restricted to some areas, or is known to concentrate heavily in some areas, then using the multiple groups version is more appropriate.We found that a minimum of 15 groups is requires to appropriately avoid the assumption of regeneration of spatial distribution.The multiple groups formulation will allow local depletion to occur and persist through the life of each cohort.
An important feature of the modeling approach presented is that it explicitly accounts for spatial fishing effort with the gravity model (Caddy 1975).Because the effort dynamics are modeled D r a f t explicitly, it is possible to use the model to investigate questions relating to active avoidance of specific areas due to high bycatch occurrence or spatial aggregation in fishing effort.In the specific example of Pacific hake, the U.S. fleet has a strong incentive to avoid areas where abundance of bycatch species is high, despite potential high abundances of the target species.This effect can be modeled by linking the E pot,r to the abundance and distribution of bycatch species.Spatial aggregation of fishing effort is also important because if fishers aggregate in areas of high abundance, it is likely that strong cohorts will be targeted disproportionately causing variation in selectivity across time.Time varying selectivity can be confounded with fishing mortality estimates leading to biases in the estimation of reference points and management targets (Punt et al. 2015;Martell and Stewart 2014).
We found that it is possible to estimate the driving movement parameters of the model using spatial catch at age information collected throughout the migration cycle of the species being modeled.Data were generated and estimated using both single and multiple groups versions of the model.This procedure resulted in unbiased parameter estimates for the single group simulation scenarios and up to 30% median relative error in parameter estimates for multiple group scenarios.
It is unclear what are the causes for biased parameter estimates in the multiple groups version.
However, it does seem that the multiple groups version of the model has higher data requirements, which is understandable given the higher degree of complexity in the model.We have not tested the performance of this model assuming catch at age from more than five fishing grounds, but it is possible that less aggregated data would result in better resolution of the model parameters.However it is also important to consider that less aggregated data will probably have higher observation error levels and that might also impair the estimation of the model parameters.
It is possible that the relative error estimates obtained in the simulation-estimation analysis are over optimistic because the same model structure was used for simulating data and for estimating parameters.This similarity is likely to have improved the realized model fit.In addition, the pa- .However, when we attempted to jointly estimate the movement and population dynamics parameters, i.e. use the current model as a stock assessment model, there was a great amount of confounding between the estimates of productivity, recruitment deviations and the movement parameters (results not shown).Therefore, it is unlikely that the movement dynamics presented here could be integrated into stock assessment models.
A promising application of the Lagrangian model described here is its potential to be used as an operating model in closed loop simulations.Such simulations can be used to evaluate the effects of management strategies for exploited fish populations (Giske et al. 2001;Sainsbury et al. 2000).
The model can be used to represent the complex population dynamics of migratory species, as well as the variability in distribution of stocks due to intrinsic (e.g.growth, maturity) and extrinsic (e.g.environmental forcing, fishing effort) forces to the population.One advantage is that the mechanics of this model is different from that usually implemented in stock assessment models (e.g., Methot and Wetzel 2013;Fournier et al. 1998), which would yield significant differences between operating and estimation models in closed loop simulations.Similarities between operating and estimation models lead to improved performance of estimation models, which in turn result in overly favorable performance of management strategies (McClure et al. 18-21, February 2014).In reality it is unlikely that the estimation models capture all the processes that occur in nature.
A couple model extensions were not included in this paper, but could easily be added to the model dynamics.These extensions are the addition of multidimensional movement and the use of other mathematical functions to represent migration hypotheses.The model illustration presented in this study only describes movement in a unidirectional basis, that is from spawning to feeding grounds.This simplification of movement trajectory stems partly from the knowledge of the species chosen as a case study.A unidimensional model is commonly used to describe the movement dynamics of Pacific hake and not much is known regarding the population trajectory and habitat use in other dimensions (Ressler et al. 2007(Ressler et al. (2008)).However for many species, migration routes are more complex and involve simultaneous migration between spawning and feeding grounds, shallow and deep water and between inshore and offshore grounds (e.g., Misund et al. 1998;Barbaro et al. 2009;Merten et al. 2016).The current approach could be ,easily extended to D r a f t a multidimensional approach.The addition of new dimensions, however, would require the development of mathematical functions that describe movement in each dimension.The sine function presented here is a good candidate to describe cyclic movements, but linear, logistic or knife-edge functions could also be used if movement in a given dimension is thought to be permanent, as would be true for migration between nursery and rearing grounds.
We believe that the model presented here is a useful approach to model movement of migratory fish species.We anticipate that the model can be used to examine the plausibility of different movement hypotheses and to explore the possible links between fish migration and ecosystem interactions.In addition, we suggest that the model is a good candidate to be used as an operating model in closed loop simulations, especially when there is interest in evaluating the implication of migratory movement on management outcomes.

Acknowledgments
This work was funded by the NSERC Canadian Fisheries Research Network.We would like to thank Aaron Greenberg for letting us use his personal computer to perform the simulation exercise shown in this study.We also would like to thank two anonymous reviewers and Dr. Daniel Goethel for his comments and extensive review of this manuscript.constant for all years and equal to (1, 1, 1, 0.2, 0.2) for fishing grounds from 1 to 5 yearly effort scaler -fishing grounds 1-2 and 4-5 were combined when only 3 fishing grounds were considered E t,k (0.0 0.0 0.0 0.0 0.5 1.0 1.0 1.0 0.5 0.1 0.0 0.0) for fishing grounds 1-3 and (0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.5 0.3 0.0 0.0) for fishing grounds 4-5 monthly effort scaler -fishing grounds 1-2 and 4-5 were combined when only 3 fishing grounds were considered q 1 effort scaler D r a f t or continuous advection-diffusion models(Sibert et al. 1999).Box models are simpler but require the predetermination of, usually large, spatial areas from which flow (i.e.movement) is measured.The determination of such areas can be challenging, especially due to the assumption that fishing mortality is homogeneously distributed within each area and that flow in between boxes is mainly caused by migration or diffusive movement(Walters and Martell 2004;Carruthers et al. 2011).Frequently, these large spatial areas are determined by political or management boundaries or historical division of data, and do not correspond to an ecologically relevant partition of the habitat of the species being studied.Such artificial partitions can lead to violation of the box model assumptions.The Lagrangian model presented rameter estimates obtained in this study are dependent on assumptions regarding the fixed model parameters, i.e. the fisheries and population dynamics parameters.The true values of such parameters are usually not well known and estimates can change dramatically over time(Brooks et al.

Figure 2 -Figure 3 -Figure 5 -Figure 6 -(
Figure 2 -Map of the study area.The dashed lines indicate the division between fishing grounds. All three sections are structured by age, time, space (i.e. modeled area, fishing grounds and territories), and group.All model indexes, i.e. variables used to designate the model dimensions, are presented in Table 1.The age dimension, denoted as a, aggregates individuals by cohort, and spans from age one to a plus group, which encompasses individuals of age A and older.Individuals age is set to increase only when t = t 0 .
Melnychuk et al. 2017)and the collection of such data is already part of the management programs for such species, i.e. used in stock assessments.Therefore, the Lagrangian model presented can be used as an alternative tool to formalize and test the movement hypotheses based on data, even when tagging or movement track data are not available.It is important to note, sented here relies mainly on seasonal catch and age composition data to estimate its movement parameters.Catch and age composition data are conventionally available for many temperate ex-

pubs Canadian Journal of Fisheries and Aquatic Sciences
Melnychuk, M.C., Peterson, E., Elliott, M., and Hilborn, R. 2017.Fisheries management impacts on target species status.Proceedings of the National Academy of Sciences 114(1): 178-183.Methot, R.D. and Dorn, M.W. 1995.Biology and fisheries of North Pacific Hake (M.productus).In Hake: Biology, Fisheries and Markets, Springer, [S.l.], pp.389-414.Methot, R.D. and Wetzel, C.R. 2013.Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management.Fisheries Research 142: 86-99.doi:10.1016/j.Pine, W.E., Pollock, K.H., Hightower, J.E., Kwak, T.J., and Rice, J.A. 2003.A Review of Tagging Methods for Estimating Fish Population Size and Components of Mortality.Fisheries 28(10): Walters, C.J. and Martell, S.J.D. 2004.Fisheries Ecology and Management.Princeton University Press.Waterhouse, L., Sampson, D.B., Maunder, M., and Semmens, B.X. 2014.Using areas-as-fleets selectivity to model spatial fishing: Asymptotic curves are unlikely under equilibrium conditions.Diagram illustrating the sine function used to describe the cyclic movement dynamics used in the Lagrangian model.The position (y-axis) of an individual changes in cyclic waves as time (x-axis) progresses.In this diagram, we assume that a cycle lasts 12 time-steps and that movement occurs between values of 30 and 60 along the spatial scale. . . . . . . . . . . . . . . . . . . . . . . . . . ..24 2 Map of the study area.The dashed lines indicate the division between fishing grounds. . . . . . . . . 25 3 Illustration of the distribution of fish of age 5, in the month of July for the single and multiple groups versions in the absence of fishing.The multiple groups version output is shown for all groups individually and combined.The vertical dashed line represents the border between territories. . . . . . . ..26 4 Monthly representation of differences between spatial distribution of biomass for single and multiple groups versions.Vertical bars indicate relative fishing effort.Continuous and dashed lines represent vulnerable biomass for fish at age 1 and 5 throughout one year migration cycle, one panel for month).Different shades of gray indicate multiple and single group versions.The vertical dashed line represents the border between territories.Both fish biomass and effort were re-scaled from 0 to 1 for plotting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..27 5 % Relative error for estimation of key parameters of the movement model.Scenarios are indicate in the upper right corner of each plot and the scenarios composition is indicated in the plot titles.SG = single group version , MG= multiple groups version and FG = fishing grounds. . . . . . . . . . . ..28 6 Median biomass distribution for simulation and estimation models for each month of the year for scenario 11 -multiple groups, 5 areas, τ=1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..29 Diagram illustrating the sine function used to describe the cyclic movement dynamics used in the Lagrangian model.The position (y-axis) of an individual changes in cyclic waves as time (x-axis) progresses.In this diagram, we assume that a cycle lasts 12 time-steps and that movement occurs between values of 30 and 60 along the spatial scale.

Table 6 -
Pacific hake model dimensions and parameter values

Table 7 -
Simulation-estimation scenarios Scenario number Model Version catch at age error Fishing grounds λ value