FISH POPULATION DYNAMICS MODEL COMBINED WITH A STATISTICAL PARAMETER ESTIMATION PROCEDURE: APPLICATION TO THE INDIAN OCEAN SKIPJACK TUNA FISHERY

We develop an advection-diffusion size-structured fish population dynamics model and apply it to simulate the skipjack tuna population in the Indian Ocean. The model is fully spatialized, and movements are parameterized with oceanographical and biological data; thus it naturally reacts to environment changes. We first formulate an initial-boundary value problem and prove existence of a unique positive solution. We then discuss the numerical scheme chosen for the integration of the simulation model. In a second step we address the parameter estimation problem for such a model. With the help of automatic differentiation, we derive the adjoint code which is used to compute the exact gradient of a Bayesian cost function measuring the distance between the outputs of the model and catch and length frequency data. A sensitivity analysis shows that not all parameters can be estimated from the data. Finally twin experiments in which pertubated parameters are recovered from simulated data are successfully conducted.

1. Introduction.Fish population dynamics models together with parameter estimation techniques are essential to provide assessment of the fish abundance and fishery exploitation level.Their use forms the basis of scientific advice for fisheries managements.This is particularly true for tuna fisheries, which are among the most valuable in the world and subject to increasing fishing pressure and to the effects of climate changes.
Discrete age-structured models with crude representations of space are most of the time used for fisheries stock assessments [1,2].The classical data used in fishery science to calibrate models are fishing effort, catch and length frequency data.
Length frequency data are not straightforward to use.Fish of the same age can exhibit very different sizes depending of their history [3,4].Therefore, to compare the outputs of age-structured models with length frequency data, a Gaussian size distribution is generally added to each age class.However, because of non-uniform mortality over sizes, bias on growth and mortality estimates may result from this procedure [5].
Another point concerning tuna fisheries is that they are highly heterogeneous in space and time.This has a significant effect on their functioning.Important migrations of fish occur at various scales, so that fish movements have to be explicitly represented using spatialized models [6].
These are some of the main problems of current stock assessment models.It is necessary to carry on the modelling effort by proposing and testing more complex models dealing more accurately with size distributions and spatial heterogeneity.This paper follows this direction and its purpose is twofold.
First we describe in section 2 a model of population dynamics in which both size and space are taken as structure variables to account for growth, movements of fish, environmental variability and variable distribution of fishing effort.The model consists of an advection-diffusion-reaction equation.Spatial advection-diffusion models have a long history in ecology [7,8,9], but their use in fishery science has grown recently, particularly for tuna population modeling purposes [10,11,6].To model tropical tuna population in the Indian Ocean realistically, our model needs to reflect the heterogeneous distribution and movements of the population linked to the environment and fishing effort heterogeneity.Thus in the model fish movements depend on oceanographical and biological data through a habitat suitability index.Recruitment, that is to say the input of young fishes in the model, is modeled as a source term involving a nonlocal nonlinearity.The structure of the model enables a direct and simultaneous comparison with the two main types of data available for tuna fisheries: catches and size frequencies.
We assess the mathematical well-posedness of the model in section 3. We formulate an initial-boundary value problem, introduce a variational formulation and show existence of a unique weak solution.As often with nonlinear problems the proof uses a fixed-point argument.We also show the positivity of the solution.
Our second goal is to develop and test a data assimilation procedure to estimate the parameters of the model in a realistic skipjack tuna fishing simulation.Indeed one of the main objectives of tuna population modeling is to provide robust evaluations of stocks which are hardly possible nowadays for tuna fisheries in the Indian Ocean because of the lack of robust estimations of many biological parameters, such as natural and fishing mortality rates or recruitment parameters.Section 4 deals with the numerical implementation of the simulation model.Then in Section 5 we describe the data assimilation method developed for parameter estimation as well as the Bayesian likelihood approach used to formulate a cost function measuring the distance between the outputs of the model and the data.The paper ends with some numerical experiments conducted in section 6 to validate the algorithm in the case of the Indian Ocean skipjack fishery.
2. The model.The dynamics of the population of fish is described through a density function p(x, y, s, t), where position (x, y) ∈ Ω the bounded domain representing the ocean, size or length s ∈ (S 0 , S 1 ) and time t ∈ (0, T ).The number of fish of size between s 1 and s 2 at time t and position (x, y) is given by the integral s2 s1 p(x, y, s, t)ds.

ADVECTION-DIFFUSION-REACTION SIZE-STRUCTURED FISH POPULATION MODEL 3
The population density follows an advection-diffusion process in space.Let D(x, y, s, t) = diag(D u (x, y, s, t), D v (x, y, s, t)) be the space diffusion matrix and V (x, y, s, t) = (u(x, y, s, t), v(x, y, s, t)) T be the velocity field.The population density also follows an advection-diffusion process in the size variable (see section 2.1 for more details).Let d(s) denote the dispersion coefficient in size and γ(s) be the growth rate.Finally let m(s) and F (x, y, s, t) denote the natural and fishing mortality rates, and R(x, y, s, t, p) the recruitment source term (see sections 2.4 and 2.3).The density function p follows the balance law, where ∇ and div are the usual differential operators on Ω.This equation has to be completed with initial conditions and boundary conditions and ∇p(x, y, s, t) • n(x, y) = 0, in ∂Ω, ∀(s, t) ∈ (S 0 , S 1 ) × (0, T ), where n(x, y) is the unit normal vector pointing outside Ω.Homogeneous Neumann boundary conditions at s = S 0 and s = S 1 express the fact that the size of individuals can not reach values lower than S 0 or larger than S 1 .
The parameterizations of the processes involved in the time evolution of the population are described in detail in the following subsections.
2.1.Movements: Advection-diffusion in space.Diffusion and velocity in space have a physical and a biological component.The biological components depend on a habitat suitability index function, hsi(x, y, t), and its first space derivatives.The index hsi depends on temperature, T (x, y, t), and forage, F ood(x, y, t) which are input data for the model.The biotic affinities for these environmental factors are defined as and All parameters are given in Table 1, and Fig. 2 shows plots of f T and f F ood .In its most general form the index is defined as The velocity field is computed as where V phy = (u phy , v phy ) T represents the physical velocity (computed by a hydrodynamical model) and V hsi = (u hsi , v hsi ) T represents biological velocity defined by ), V hsi is proportional to length (large fish can swim faster than small ones) and to (1−hsi)∇hsi (the model transports the population towards the most suitable places for fish living according to the temporal habitat index evolution).The diffusion coefficients are defined as follows: The interpretation of such a parameterization is similar to the one given for V hsi and all parameters are given in Table 2.  12) and Table 3) is included to account for individuals having the same age but different sizes.Indeed, in a fish population individuals of the same age can often differ markedly in size [4].This variability in growth can result from many mechanisms, including genetic or behavorial traits that confer different performances to individuals, and factors such as environmental heterogeneity and variability [3].In fishery science, this variability is usually taken into account in age-structured models using a lengthat-age relation perturbed by a Gaussian noise (see, for example, [12]).The model discussed here is size-structured and uses a diffusion term in the size variable with dispersion rate d(s) to account for individuals having the same age but different sizes [13,5].The advection-diffusion term in size can be seen as the limit of a random walk model in which each individual grows with an average velocity but has at each time step a small binomial probability to grow faster or slower than this average (see [8] for more details).We consider that fish growth follows a Von Bertalanfy curve: where s mat is the minimum size at maturity, the fecundity rate f r(s) is given by and the weight w(s) of a fish of size s by We use a Beverton and Holt [3] stock-recruitment relation and obtain, where l 1 (S0,sr) is the usual characteristics function and s r is the maximum size of recruitment.

2.4.
Natural and fishing mortality.The mortality rate is split into size-dependent natural mortality m(s) [14] and a fishing mortality rate F (x, y, s, t).The fishing mortality rate is defined as the sum of the N f fishing mortality induced by each fleet, The mortality rate induced by each fleet is described by the following equation: where q f (s) is the size-dependent catchability coefficient for fleet f (that is the probability for a fish of size s to be caught by a unit of fishing effort of fleet f ), and E f (x, y, t) is the observed fishing effort.
3. Mathematical well-posedness.In this section we prove existence and uniqueness of a positive weak solution to the model.

Functional spaces.
Let us introduce some functional spaces.The study is conducted on the open set Q = Ω × (S 0 , S 1 ).Let T < ∞ be a fixed time.
The population density, p, is considered as an element of the functional space H = L 2 (Q), whose Hilbert space machinery is convenient to use.H is equipped with the scalar product pqdsdxdy, and we denote by ||.|| H the induced norm.We also consider the separable Hilbert space defined by H 1 = H 1 (Q) and equipped with the scalar product (pq + ∇p.∇q + ∂ s p∂ s q)dsdxdy.
We denote by ||.|| H 1 the induced norm on H 1 .We will also have to consider the Banch space and L ∞ (0, T, H) is the space of functions L ∞ in time with values in H, equipped with the norm, Similarly C([0, T ], H) is the space of continuous functions on [0, T ] with values in H. Further, C([0, T ], H), L ∞ (0, T, H) and L 2 (0, T, H) are Banach spaces.Classically H ′ denotes the dual of H and (H 1 ) ′ the dual of H 1 .When H is identified with its dual, we have the scheme where each space is dense in the following and the imbeddings are continuous.Let us denote by W (H 1 ) the Hilbert space ) is a.e equal to a continuous function from [0, T ] to H. Moreover we have the following continuous imbedding, Proof.See, for example, Dautray and Lions [15] 3.2.Assumptions on the data and preliminary transformation of the system.The mortality rates are assumed to satisfy If we assume that the input temperature and forage fields, T (x, y, t) and F ood(x, y, t), are positive and regular enough, it appears clearly from section 2 that • u(x, y, s, t), v(x, y, s, t) are differentiable with respect to x and y, respectively and It is also clear from section 2 that We also assume that the initial distribution p 0 (x, y, s) satisfies To prove our existence-uniqueness result, it is convenient to perform a change of unknown function: p satisfies (1)-( 4) if and only if p = e −λt p is a solution to the same equations where −(m + F )p is replaced with −(m + F + λ)p in Eq. ( 1) and the recruitment R(x, y, s, t, p) is replaced by B(x, y, t, p) = S1 smat f r(s)w(s)p(x, y, s, t)ds.
In the remaining part of the mathematical analysis, this change of unknown is implicitly done and we omit the p notation.The constant λ will be fixed to a convenient value below.Moreover, the possible nullification of the term k B e −λt + B(x, y, t, p) invites us to define This formulation is being used in the following.We will show that if the initial distribution, p 0 is nonnegative then p ≥ 0 a.e. in Q × (0, T ); thus the two formulations are equivalent.

3.3.1.
The bilinear form a(t, p, q).Formally multiplying Eq. ( 1) by a function q and integrating by parts on Q leads to the definition of the following bilinear form.For p, q ∈ H 1 let us define Lemma 3.2.For a.e.t ∈ (0, T ), a(t, p, q) is continuous on H 1 × H 1 , and for λ large enough, a(t, p, q) is coercive on H 1 .There exist two constants Proof.The proof is classical and we omit it.

3.3.2.
The nonlinear operator R. In this section we show that the recruitment term R(x, y, s, t, , p) (cf.Eqs. ( 13)-( 16)) allows us to define a Lipschitz continuous operator R on L 2 (0, T, H).

Weak solutions.
We can now give the definition of a weak solution to system (1)-( 4).
Then we can prove the following result.
Existence and uniqueness.The proof consists mainly in defining a nonlinear operator Θ by freezing the nonlinear term Rp and applying Banach fixed-point theorem to Θ.The fixed point is the desired solution.
Step 1.Let p be fixed in W (H 1 ) and in Eq. (30) let us replace (Rp, q) H by (Rp, q) H .The problem becomes linear in p and admits a unique solution (e.g, [15]).This solution defines an operator Θ on W (H 1 ), Θp = p.
Step 2. Let us show that for T sufficiently small Θ satisfies the following properties: 1. Θ leaves invariant the ball, Taking q = p as test function in (30), integrating on [0, t], using the coerciveness of a and Cauchy-Schwarz inequality, we obtain For all α > 0, Young inequality leads to and choosing α = 2C 2 gives Then using Eq. ( 28) we obtain 2. Θ is a strict contraction on B r , there exists 0 Let p 1 = Θp 1 and p 2 = Θp 2 .Substracting the two associated Eq. (30), taking p 1 − p 2 as test function and again using the coerciveness of a, Cauchy-Schwarz and Young inequality leads to Since p 1 (0) = p 2 (0) = p 0 , we deduce and Then for Λ 2 T 2C2 < 1, Θ is a strict contraction.
Step 3.For T small enough, by Banach-fixed point theorem Θ admits a unique fixed point which is the desired solution on (0, T ).Since T does not depend on p 0 , the same procedure can be applied on (T, 2T ), ... until a solution is found on the desired time interval.
Positivity.Let p 1 ≥ 0 be given in W (H 1 ), and let us define the sequence (p n ) n≥1 by Θp n = p n+1 .Let us prove that p 2 is non-negative: Taking p − 2 = max(0, −p 2 ) as test function in (30) leads to and therefore to 1 2 0, and p 2 ≥ 0. An induction then shows that p n ≥ 0, ∀n ≥ 1, and since the sequence converges to the solution p, this latter is non-negative.

4.
Numerical treatment of the model.In the approximation procedure of the model a centered finite difference discretization is used.Equation ( 1) is solved on a grid with a spatial resolution of 2 degrees, i.e. ∆y = 120 nautical miles in the latitudinal direction and ∆x = ∆y cos(θ) in the longitudinal direction (θ is the latitude angle), a discrete length step, ∆s of 4 cm, and a discrete time step ∆t of one day is used.The discretization points are denoted by (x i , y j , s l , t n ) with i ∈ [1 : I], j ∈ [1 : J] (assuming here for simplicity that the domain Ω is rectangular), l ∈ [1 : L] and n ∈ [1 : N ].In what follows, p n i,j,l denotes the numerical approximation of p(x i , y j , s l , t n ).
Several difficulties arise in the computation of the solution to Eq. (1).First the numerical scheme has to be very stable because of possible strong variations in space and time of the advection and diffusion coefficients, u, v, D u , D v .Moreover the numerical solution of Eq. ( 1) is to be used in a numerical function minimization procedure to obtain estimates of model parameters.Therefore the solution algorithm must be fast because the model and its adjoint may have to be solved hundreds of time.Moreover, the function minimization algorithm may test parameter values that do not necessarily guarantee numerical stability.
The selected scheme combines a splitting method [16,17] and the use of the MUSCL scheme for advection terms (monotonic upstream centered scheme for conservation laws [18]).At each time step, given an approximation p n of p(x, y, s, t n ), the computation of p n+1 from p n is achieved through four steps.The advectiondiffusion equation in the x variable is integrated first, on [t n , t n+1 ]: It results in a first approximation p n+1,1 .Then the advection-diffusion equation in the y variable is integrated on [t n , t n+1 ] starting from p n+1,1 : It results in a second approximation p n+1,2 .Then the advection-diffusion term in the s variable is integrated on [t n , t n+1 ] starting from p n+1,2 : It results in a third approximation p n+1,3 .Finally mortality and recruitment are integrated on [t n , t n+1 ] starting from p n+1,3 : It results in the final value p n+1 .In each of the first three steps, diffusion is treated implicitly in time, and the MUSCL scheme is used for the advection term.For example, the discretization used to solve Eq. (31) can be written as follows: where and otherwise, 5. Parameter estimation.In this section we describe the data assimilation algorithm developed to estimate the parameters of the model.

5.1.
The outputs of the model corresponding to the data.Total catches in weight as well as length frequencies of the catches are computed and compared to observations to estimate the parameters of the model.In each cell (i, j) of the grid, where during month m (30 days), the fishing effort is nonzero, catches of fleet f are computed as follows: and length frequencies as We use the maximum of posterior distribution method [19], which involves minimizing the sum of the negative log-likelihood of the data plus the log of prior density functions.We assume that the observation errors for catch data follow a log-normal distribution.Therefore the contribution of total catches to the negative log-likelihood hal-00496619, version 1 -1 Jul 2010 The observation errors for length frequency data are assumed to be normal and the contribution of frequency data to the negative log-likelihood reads The negative log of prior density functions for the parameters is where K 0 n are the reference a priori parameters given in Tables 1-4.The cost function to be minimized is the sum of those three terms: The parameters have different units and orders of magnitude.To avoid any numerical difficulties that might arise from this during the minimization, we adimensionalize the parameter vector K, dividing each parameter K i by its first guess a priori value K 0 i .Let D = diag(K 0 i ) then the adimensionalized control vector is k = D −1 K.Such an adimensionalization procedure can be regarded as a preconditioning for the minimization.The final cost function is and the a priori reference adimentionalized parameter vector is k 0 = 1.

Optimization:
Computing the gradient with the adjoint model.To minimize the cost function j, we used the quasi-Newton algorithm implemented in the n1qn3 Fortran subroutine of Gilbert and Lemaréchal [20].The computation of the gradient of j with respect to control variables is required at each step of the minimization.This gradient results in one integration of the adjoint model.The adjoint code was obtained using the automatic differentiation program Odyssée [21,22], which is an efficient tool for deriving adjoint codes since it enables the automatic production of adjoint instructions.However, codes produced by automatic differentiation do not usually use computer memory in a very efficient way.Saving the direct model trajectory is the major problem.A differentiation program has to follow systematic methods to provide the evaluation trajectory.Thus Odyssée systematically uses a local calculation and storage technique for the trajectory.Automatically differentiating a 3D model and using the adjoint code directly seems impossible for the moment.Thus the code generated by Odyssée had to be improved manually.A Taylor test was then conducted to compare the exact derivatives computed by the adjoint code to a finite difference approximation.Generally speaking, one aims at verifying that for any direction of perturbation δk.We present in Table 5 the result of such a test.As ǫ becomes smaller, one observes that the ratio r(ǫ) first tends linearly towards 1 up to ǫ = 10 −6 , which is the optimal value for a finite difference computation.Afterwards, the substraction of close floating-point numbers leads to a large cancellation error, which dominates the truncation error coming from the computation of Table 5. Result of a Taylor test ǫ r(ǫ) 10 −1 1.010756398 10 −2 1.001203422 10 −3 1.000115118 10 −4 1.000015861 10 −5 1.000001463 10 −6 1.000000518 10 −7 0.999993143 10 −8 0.999985731 10 −9 0.999847627 10 −10 0.997712868 the gradient by the finite difference method.A Taylor test with such a numerical behavior of the ratio r(ǫ) is said to be correct.It verifies that the adjoint code provides an exact computation of the gradient.Inititial conditions are chosen to be homogeneous over the space grid with a size distribution p 0 (x, y, s) = 0.1e −0.5s .
(47) This distribution assumes that the population is dominated by small organisms.Using these initial conditions a spin-up run of 6 years is conducted in order to reach an experimental and numerical fixed-point where mortality processes balance the recruitment process and the total biomass slowly varies around a mean value during the year.The final distribution at the end of the spin up period provides the initial distribution for the standard run.
Since skipjack tunas inhabit the surface layer of the ocean, the model is forced with monthly velocity of oceanic surface currents, sea surface temperature and forage fields (Fig. 2).Velocity and temperature fields are outputs of the ocean general circulation model OPA,1 whereas forage fields are outputs of a size structured model representing the energy flow in marine ecosystems from zooplankton to organisms of the size of tuna forage [23].Frequency distributions over the whole grid and the whole year of temperature and forage values are also plotted.
Figure 3 shows a contour plot of the index hsi for the month of march.The strong south-north gradient of the index in the lower half of the map is typical of the area.Low temperatures are not suitable for tuna, which stay in the upper half of the map as shown on Fig. 4.
To test the possibility of estimating some parameters of the model from standard fishing data, we conduct in the following section numerical experiments with a synthetic data set computed by one simulation of the model.All parameters are set to their reference a priori values.Moreover for the sake of simplicity the two mortality parameters m and q are assumed to be size-independent; that is constants with values m = 4.2438 10 −8 s −1 [6] and q = 6.43 10 −8 s −1 .Only one fleet is considered (N f = 1), and the fishing effort is assumed to be constant and homogeneous during the year on an area roughly corresponding to the real fishing areas of the purse seine fleet in the Indian Ocean (see Fig. 5).With this configuration, a data set is computed following Eqs.(39) and (40).6.2.Sensitivity analysis.We conduct a sensitivity analysis in order to identify the most important input parameters whose changes impact the most the outputs of the model (catches and length frequencies).As a measure of the outputs of the model we consider the quantities Then the vector of relative sensitivities of these quantities to variations of the input parameters computed at the point k are and Each of these relative sensitivity vector requires the computation of a gradient which is easily obtained with one integration of the adjoint model.Table 6 shows the relative sensitivity vectors, s C (k 0 ) and s Q (k 0 ), computed with the initial a priori parameter vector k 0 .Globally the different relative sensitivities are quite low, indicating that it may be difficult to estimate correctly all the parameters of the model using the two types of fishery data which are generally available.In the remaining part of this paper we will not try to estimate parameters which have the lowest sensitivities (10 −3 , 10 −2 ).These parameters are • D min , D max related to diffusion in space, • d 1 , d 2 , γ 1 and γ 2 related to growth, • a f and s mat related to recruitment.Interestingly the relative sensitivities s Q corresponding to variations in b f , k B and a w (see Table 6) are exactly equal up to a sign.This comes from the formulation of recruitment in the model, which from Eqs. ( 13)-( 16) can be rewritten with obvious notations as An inconsistency in the formulation of the inverse parameter estimation problem appears clearly.The 3 parameters, b f , k B and a w can not be determined independently, since for example an increase in k B can also be interpreted as a decrease in b f or in a w .For this reason in our identification experiments we keep b f and a w fixed to their reference values and only try to estimate k B .Moreover since the length/weight parameters a w and b w are well known we also do not select b w for the parameter estimation formulation.
Although the 2 mortality parameters m and q do not correspond to very high sensitivities, we will try to estimate them since they really are badly known.
Finally the chosen formulation includes 11 parameters to be estimated: • movements parameters : k hsi , p T , p F ood , α T , T 0 , K F ood , u hsi0 • recruitment parameters : b 0 , k B • mortality parameters : m, q The assumed variances are as follows: σ C = 0.1 and σ Q = 0.01.This provides a good balance between the two terms of j.In the experiments conducted, the convergence criterion is where ǫ is a small value fixed to 10 −5 .

Experiment 1.
A first numerical experiment was conducted to assess the capacity of the parameter estimation algorithm to distinguish between low (or high) recruitment and high (or low) mortality rates on the one hand and natural and fishing mortality rates on the other hand.Therefore in this optimization only the four parameters b 0 , k B , m, and q can vary, the others being fixed to their reference a priori value used to simulate the data.Different first guesses for the parameter vector were obtained by perturbing these four parameters within reasonable range (up to 50% of their reference value).All the corresponding optimizations converged

Conclusion.
We developed an advection-diffusion size-structured fish population dynamics model and applied it to simulate the skipjack tuna population in the Indian Ocean.The model is fully spatialized and movements are parameterized with oceanographical and biological data.Thus the model naturally reacts to environmental and climatic changes.We have formulated an initial-boundary value problem and proved its mathematical well-posedness.We then discussed the numerical scheme chosen for the integration of the simulation model.From a modeling point of view, this study, is to our knowledge, the first one in which space and size structure of the population are fully taken into account and in which both mathematical and numerical difficulties were dealt with in a rigorous manner.
In a second step we addressed the parameter estimation problem for such a model.With the help of automatic differentiation we derived the adjoint code which enabled us to compute the exact gradient of a Bayesian cost function measuring the distance between the outputs of the model and catch and length frequency data.Thanks to the size structure of the modeled population the outputs of the model can be naturally compared to length frequency data.A sensitivity analysis showed that not all parameters could be estimated from the data.Finally twin experiments in which pertubated parameters were recovered from simulated data were successfully conducted.This point is particularly crucial since one limitation of the model lies in the choice to be made for different parameters value, or even in the choice to be made in the type of functions of temperature or forage parameterizing the habitat.The numerical experiments conducted demonstrate that fishing data can be used to estimate these parameters accurately.
This study is an important first step towards the assimilation of real observed fishing data in the model which is under progress.The mathematical and numerical tools which have been developed and validated will be extended to confront the model with tagging data which should bring more information and enable the estimation of several supplementary parameters such as growth and movements parameters.Developing a tool using tagging data is indeed especially timely, since no reliable stock assessment can be conducted at present for the skipjack tuna in the Indian Ocean and since a large-scale tuna tagging program in the Indian Ocean (IOTTP) has recently started and an important tag-recapture data set will be available in the coming months.

) 5 . 2 .
The cost function.The parameters of the model are denoted in what follows by K ∈ R Np where N p is the number of parameters.K is being estimated in a Bayesian context by computing the mode of the posterior density function of the parameters knowing the data.

6 .
Numerical results.6.1.The simulation set-up.The standard run consists in a one-year simulation for the Indian Ocean.The spatial numerical grid used is shown on Fig. 1.Sizes of simulated skipjack tunas range from S 0 = 0.4 m to S 1 = 1.2 m.

Figure 1 .
Figure 1.The 46 × 32 numerical grid used to integrate the population dynamics model on the Indian Ocean (earth in black, ocean in white).

re s p o n s e fu n c ti o n v a l u eFigure 2 .
Figure 2. Response functions, f T and f F included in the hsi formulation corresponding to temperature (left) and forage (right).Frequency distributions over the whole grid and the whole year of temperature and forage values are also plotted.

Figure 3 .
Figure 3. Contour plot of the index hsi for the month of March.

Figure 4 .
Figure 4. Contour plot of the population density function summed over size classes at the end of the month of March.

Figure 5 .
Figure 5. Distribution of the fishing efforts used in the simulations.F = 0 in the gray area and F = 1 in the black area.

Table 1 .
Habitat suitability index parameters

Table 3 .
Parameters for growth and dispersion in size

Table 6 .
Relative sensitivities of catch and length frequency data Identification experiments.An essential validation step to perform before assimilation of real observed data is to conduct twin experiments.Synthetic data are produced by the model using the first guess parameter vector k 0 .To fully test the possibility of recovering the selected parameters from the synthetic data, no penalty term is added and the cost function reduces to j