The purpose of the model is to predict spatial patterns of abundance for the spawning migration along the Norwegian coast. This model provides discrete estimates across the model area that can be used to compare against concurrent observations. Furthermore, this model is designed to improve estimates when observations become available. A brief schedule of the main operations is presented in Table 1. The movement of individuals is modelled by changes in orientation and horizontal speed. In particular, the response of individuals to the Norwegian Coastal Current (NCC), along with temperature and depth gradients are modelled. The individuals also utilize knowledge of previous states, such as orientation angle, when moving.

Estimates of environmental conditions were loaded from SINMOD, a physical oceanographic model that is based on the primitive Navier–Stokes equations, and uses a z-coordinate grid (Slagstad and McClimans, 2005). The configuration used has 970 × 635 horizontal grid cells with 4km resolution and centres on the Norwegian Sea. The model is divided into 34 vertical depth layers. The IBM developed in this paper modelled individuals in a 2D environment where position was updated on a continuous horizontal plane in a Lagrangian approach (Figure 1). Environmental variables were calculated based on assumptions of the herring's depth preferences, described in section 2.4.1. State variables were updated at discrete time steps of Δt = 4h. Temperatures and current speeds were extracted from SINMOD output from 2015 to 2020, and along with the bathymetry field of the model area, drove changes in fish movements.

The model was developed in MATLAB, which is a matrix-based programming language that is suited for iterative analysis involving numerous matrix operations. Below, the development of the model is outlined in regards to the model system, the main equations, parameters and state variables. Thereafter, we explore tuning and comparison against observations of the real system.

When referring to vectors that can take on continuous values, x and y indices will be used, while the index j will indicate the discrete linear index of a grid cell, ranging from 1 to the number of elements in the grid. Finally, boldface characters denote vectors.

The individual state variables used were position p, orientation angle θ, horizontal speed r_b, and horizontal speed offset r_o. Position was updated at each discrete time k, with time step Δt:

where:

where p is a vector [pxpy]T with the x and y coordinates of the individual in the continuous space of the model grid and v_b is a vector [vbxvby]T with the horizontal velocity components of an individual fish in the x and y directions, based on behavioural cues. Similarly, v_c is a vector [vcxvcy]T with the horizontal current velocity components in the x and y directions. The superscript T denotes the transpose of the vector. The spawning migration proceeds counter to the NCC (Slotte and Fiksen, 2000). This is modelled by the term −Φv_c that adds a counter-current component to the horizontal speed controlled by the parameter Φ. An individual's realized swimming velocity v_f is composed of the counter current term and behavioural responses from v_b. This formulation demands individuals respond to the prevailing current with higher priority, relative to other cues. To prevent unrealistic dynamics in the first-order approximation of velocity (v_f + v_c), the short Δt of 4h was used. The vector v_b was calculated as:

where r_b is the horizontal speed of the individual in m s⁻¹ and θ the orientation angle according to gradient cues. As indicated before, when ||v_c|| approaches zero, the individuals approach a speed of r_b with the orientation angle θ. The angle θ was updated as follows:

where α is a weighting parameter, and G is the angle of the vector inputs calculated from near-field gradients, as explored in section 2.4.2. It is likely NSSH base their movements on comparison of conditions from previous experience and present information when calculating their new orientation (Fernö et al., 1998). To account for this, α acted as a low-pass filter, avoiding erratic changes in θ, similar to a formulation by Føre et al. (2009). Furthermore, r_b was calculated as a random process with a deviation r_o from the cruising speed:

where r-b was the cruising speed in m s⁻¹ of the individual. The value was calibrated during the tuning procedure in section 2.5. The offset r_o was then calculated as a Gauss-Markov process with exponential auto-correlation. This meant r_o and r_b were correlated with recent values. The offset was included to model randomness in the horizontal speed. It was updated as follows:

from a normal distribution N with zero mean, a standard deviation σ in swimming speed and shaping parameter for auto-correlation β.

The spawning migration follows overwintering, a stationary period characterized by slow swimming speeds and low energy use. NSSH mainly spawn in southern regions of coastal Norway before migrating westwards to feed over summer. The spawning migration begins early to mid-January and spawning usually begins in late-February/early-March (Dragesund et al., 1980). Temperature and current data for the period from mid-January to the end of February (2015–2020) were loaded from SINMOD, along with the bathymetry. Below, the responses that play a role in the spawning migration are outlined.

To develop the individuals responses to environmental information, the parameters [Φ α w]^T in Equations (2), (4), and (9) were tuned and then analysed. Simulations with 2017–2020 environmental values and 4 separate initialization scenarios were investigated. In order to model reasonable responses of individuals to their environment the model was tuned using a numerical optimization algorithm that took parameter values as input and minimized the deviation between the model and observed distributions at specific time stamps. NSSH survey data from the Norwegian Institute of Marine Research was used for this purpose (Salthaug et al., 2020). The survey data combines acoustic and trawl data to estimate abundance in predefined strata areas. To ensure consistency, 2017 and 2018 estimates were used for tuning, when the dates sampled were in the last two weeks of February. 2019 and 2020 data provided an independent data set to validate the optimization.

The model output was compared with 2015–2020 catch data from Norwegian Fisheries Directorate in 5 day ranges. The catch data used were the x and y starting locations of trawling and associated catch weight in kg. The main purpose of this comparison was to investigate where and when the model deviates from observations and what this reveals about the capacity for the IBM to resemble realistic NSSH distributions. The initial model distribution was fixed in each year to focus on relative comparisons.

Observation and model densities were allocated 5 day windows, where catch and model positions and values were described by their centre of mass. These were described as latitude and longitude centre points. Then, model output was analysed using geospatial comparisons, as described in Woillez et al. (2007). The global index of collocation (GIC) indicates how geographical distinct two ditributions are. It is based on the distance between the centre points and the variance around these centres (inertia). A value of 0 means there is no overlap, whilst a value of 1 implies identical distributions. The average GIC value for each year was used to score the model and observation overlap. In addition, RMSD values indicated the average error for the year. Using GIC and RMSD indices illuminate where the model and observations disagree. It also offers insight into potential limitations of comparing model and observation output.

Of the 4 simulation scenarios run, 3 scenarios produced reasonable migration patterns. S₁, S₂, and S₃ performed relatively well (Supplementary Material). S₄ is far from continental slope, which is vital information for orientation and therefore couldn't minimize Equation (12) properly. The average Φ value, produced by the tuning process for the four scenarios, was quite high and displayed higher variability (0.78±0.2) compared to α (0.59±0.08) and w (0.65±0.1). This illustrates the sensitivity of the response to the current in relation to starting point. The relatively high α value also shows the importance of fish retaining knowledge of previous states in the migration. These simulations highlight the centrality of the continental slope as a landmark in the migration and how individuals are likely to utilize it. Further comparisons from 2015 to 2020 were made using the fitted parameter values based on S₂ values (Table 3).

The model formulation required many steps of refinement to give reasonable output. In Equations (1) and (2), the response of the individual was formulated to account for the physical response directly against the prevailing current. Removing the counter-current response led to passive drifting northwards. Running simulations with environmental variables sampled close to the surface resulted in similar drifting patterns, as the magnitude of v_c is very high close to the surface, especially around the Lofoten basin (Figure 4C). A formulation that incorporated vertical conditions was important to model, as in Equations (7) and (8).

It is difficult to decouple the low resolution effects of local environmental values on model states, as they are spatially correlated in the horizontal plane. However, the combined effect of reacting to gradients from deep bathymetry, and current patterns seemed the most consistent properties amongst simulations (Figure 4). This study shows that memory- and gradient based reactive mechanisms can be used to model the migration of pelagic fish species such as NSSH (Fernö et al., 1998).

The energetic states of individuals were not included in this study, but may provide more insight into variations in observed patterns. Further work in this project aims to incorporate Dynamic Energy Budget theory, that can be used to model energetic costs of migration (Kooijman, 2000). Thus, state-dependent selection of spawning grounds can be explored, testing the effect of body length and condition on choice of area (Slotte and Fiksen, 2000). Natural mortality may also be included here. Interactions between individuals was explored, such as attraction and repulsion, but deemed difficult to model at the 4 km² scale. There may be a case to model interactions between schools, where information transfer could potentially give access to novel information, but this is beyond the scope of the project at this stage.

It should be mentioned that the ocean models that force the behaviour of the IBM are subject to their own limitations. Ocean model outputs have uncertainty caused by limits in model resolution, our knowledge of the processes resolved by the model, uncertainty in initial values, boundary conditions, parameter values, and inaccuracies in numerical implementations (Lermusiaux et al., 2006). Additionally, climatic fluctuations can play an important role in survival of recruits and thus, biomass estimates (Toresen and Østvedt, 2000). Therefore, understanding how variability in climate can manifest as variability in individual state and parameter values is important moving forward.

The initialization procedure presented above functioned to generate a generic distribution that was applied to all model scenarios. This provided a realistic initial distribution as input for the tuning procedure. In reality, there is much uncertainty around the initial distribution of NSSH. This may be explained by changes in winter stay areas, which has historic precedent (Dragesund, 1970). Thus, in future work, we shall model initial distributions based on prescient knowledge of winter stay areas and recent catch observations preceding the model run.

Tuning the migration model against survey data allowed for high resolution, time-sensitive estimates of fish densities. There was good agreement between tuned models and independent datasets. Nevertheless, developing the means to calibrate the model was challenging. One difficulty was the progression of the migration in comparison to surveys. The survey starts in southern Norway and progresses northward, the inverse of the NSSH migration. This meant estimates at the beginning and end of survey largely consisted of acoustic absence. In the course of this work, many objective functions for tuning of different resolutions were provided for the individuals and there is room for further experimentation here.

The model was trained using data that has its own biases (e.g., gear selectivity) which may bias the model tuning as a result. Additionally, fisheries data is sparse in space and time. Immigration and recruitment will also sporadically alter fish abundance. As such, one model realization cannot encapsulate the full variation in fish distributions. Adding noise terms to model parameters can be used to produce an ensemble of estimates through Monte Carlo experiments. State and parameter estimates can be obtained effectively in this way (Ward et al., 2016).

The catch data served to judge the performance of the model across the time frame and highlight discrepancies. This shows what the model can and cannot resolve, and at what scales. The physics- based model used assumptions of responses to the prevailing current and gradients to reproduce the likely migration at large-scale. In contrast, catch data reflects which locations were targeted based on market value, quotas, weather conditions, etc. This is variability the model cannot capture. The model output is useful in predicting relevant catch locations in space and time at low resolutions (Figures 9, 10). It gives insight into where the development of model densities is matched with fishing activity and where it does not. The output will have operational use when it is combined with real-time fisheries-dependent data, as described below. The novel use of real-time observations can give time-sensitive estimates that are more reliable than use of model or observation values in isolation.

In future, information on vessel activity will be incorporated to improve observation quality. Classifying vessels in terms of search vs. fishing activity can provide fine-grain, continuous information on the behaviour of fishermen. Thus, one can improve coverage of observations and incorporate absence data. Finally, using a data assimilation procedure, we can improve estimates by combining a model, such as the one developed here, with observation data. Such estimates require expansion of the model to include uncertainties in initial conditions, environmental variables, parameters, etc. We have developed one realization of the model here, but modelling such uncertainties produces multiple realizations, from which one can calculate the most probable estimate of the model state (Evensen, 2009).