Drift Indices Confirm That Rapid Larval Displacement Is Essential for Recruitment Success in High-Latitude Oceans

Larval drift is a key process for successful fish recruitment. We used Norwegian spring-spawning herring (Clupea harengus) as model species to investigate the relationship between larval drift and recruitment. Larval drift indices were derived from simulations based on survey observations between 1993 and 2016. We show that forward simulated larval drift indices have an important positive relation to recruitment success. The relationship demonstrates elevated recruitment when larvae relocate rapidly northwards toward the Barents Sea. Negative or low larval drift indices coincide with only weak recruitment emphasizing limited survival in years with enhanced larval retention. Hence, with this work we combine drift model outcomes refined with survey data indicating that more extensive larval drift is an important component in population dynamics for high-latitude small pelagic fishes. However, larval displacement alone represents only one among many controlling factors but may offer possible predictions of the probability of higher or lower recruitment in the short term. The applicability of the drift indices is adaptable in all world oceans and all marine organisms that occupy planktonic life stages exposed to dynamic ocean currents. The study demonstrates how larval drift indices help to identify larval transport or retention to be crucial for population replenishment.


INTRODUCTION
The successful persistence of a fish population is ultimately dependent on life cycle and history closure, i.e., connectivity between each life stage, spawning adults, eggs, larvae, juveniles on nurseries and recruitment of juveniles back to the adult population (Petitgas et al., 2013). In many fish populations, adults will migrate "up stream" of suitable juvenile nursery areas so that eggs and/or larvae will be delivered to the nurseries in a suitable timeframe (Parrish et al., 1981). This window is often dictated by the development time of the early life history stages and the conditioning of a nursery area prior to the arrival of the juveniles. The drift of the eggs and/or larvae is often very influential on the success of supply of juveniles to nurseries and thereby linked to the recruitment success of the population (van der Veer et al., 2015). Life cycle closure can be characterized by the "migration triangle" (Harden-Jones, 1968), however, this framework does not have a provision for "retention of larvae" as presented by Iles and Sinclair (1982). The first allows the offspring to drift with currents away from spawning grounds toward their nurseries. The second retains the offspring in the vicinity of the spawning grounds, at least for a while, and so provide an imprinting mechanism for the return to specific spawning areas (philopatry, natal fidelity) (Sinclair, 1988). Both dispersal mechanisms are assumed to benefit larval survival depending on the physical features of the location in question (Houde, 1989).
The origin of larval drift as a key process in population dynamics of fish was firstly presented by Hjort (1926). In 1900, Johan Hjort observed early cod (Gadus morhua L.) larvae drifting in currents over hundred kilometers away from spawning grounds. With the resulting strong year-classes Hjort assumed that this larval relocation is essential for elevated recruitment. He also described that an aberrant larval drift into unsuitable nurseries may increase larval mortality rates. In contrast, instead of larval drift, retention seems to be the key when larvae instead remain at their spawning areas (Bakun, 1996). Larval retention is defined as the contrast to larval drift since larvae are kept in a certain location (Cowen et al., 2003). However, the application of the word retention is vague because it is used in several contexts. The term has been used to mean the larvae remain in one particular area (Swearer et al., 1999). In other ecosystems like important upwelling ecosystems (García-Reyes et al., 2015), larval retention is regarded as a process keeping the larvae on the shelf even though larvae can be transported alongshore (Churchill et al., 2011). Alternatively, larvae can also be retained in mesoscale eddies, even though non-stationary eddies may transport larvae into other locations (Condie and Condie, 2016).
Along the Norwegian coast, especially on the wider bank areas on the continental shelf, stationary and mobile mesoscale eddies form (Saetre, 1999). The stationary eddies are assumed to play a fundamental role in the successful recruitment of Norwegian spring-spawning (NSS) herring (Clupea harengus L.). Here larvae are suggested to be retained together with their prey on the banks (Saetre et al., 2002). This principle forms the basis of the hypothesis that retention of larval herring is dominant over the first few months of larval life (Sinclair and Power, 2015). However, there are opposing observations and hypotheses. For example, in years where early larval stages of NSS herring are exposed to a strong northward displacement this may result in strong years-classes, as was described for the period 1959period -1965period (Dragesund, 1970. In other years during this period, where larvae were retained and concentrated in specific areas, only weak year-classes occurred. The northward advection theory is further supported by recruitment success being positively linked to early spawning within a year (Husebø et al., 2009;Vikebø et al., 2010;Slotte et al., 2019) as well as to the intensity of the Norwegian Coastal Current (NCC) (Skagseth et al., 2015).
Over the last three decades, the use of biophysical modeling to identify drift patterns of marine larvae has increased more than fivefold (Swearer et al., 2019). Biophysical modeling allows particle tracking that provides information about larval dispersal which previously had to be the subject of complex and challenging field studies. Tracking NSS herring larvae in drift simulations indicated that most larvae drift was within the NCC, resulting in their arrival on their presumed nurseries in the Barents Sea plausible (Vikebø et al., 2010). However, parametrisations in drift models involve assumptions of presumed positions and times of the year a larva may have occurred in the field. To counteract this shortage, survey observations can provide a meaningful add-in to refine drift modeling outcomes (Espeland et al., 2015;Kvile et al., 2017).
The main objective of this work was to generate larval drift indices from survey observations that can be adapted to similar types of drift studies worldwide. In addition, as there is obviously some debate as to whether retention or drift results in a strong year class, further investigation is required. This clarification is especially important for oceanic, high-latitude species where the spawning is closely linked to the pulsed production of suitable prey (early copepod stages) (Sundby et al., 2016). We use NSS herring as model species and by combining information from field surveys within a biophysical model we simulated back-and forward drift of NSS herring larvae between 1993 and 2016. This period covers some successful recruitment years, but especially after 2004, where a strong environmental change was observed, there were 10 years where recruitment was particularly low . From these simulations, larval drift indices were estimated and related to a recruitment index used in a Bayesian framework to link drift strength and direction with differing levels of recruitment (Churchill et al., 2011). We address three main questions: Which dispersal pattern dominates during the drift phase of NSS herring larvae in the years between 1993 and 2016? Which dispersal patterns can be linked to the different levels of recruitment? Do and how do the larval drift indices relate with NSS herring recruitment success?

MATERIALS AND METHODS
To answer these questions field-based data were combined within advanced drift simulations. The simulations utilized a hydrodynamic dataset which allowed both back-and forward drift simulation to identify larval drift patterns for the years 1993-2016. Every effort was made to base the study on recommended practices for modeling physical-biological interactions (North et al., 2009). Any constraints in the setup of the drift simulations are further addressed in the discussion.

Using Existing Knowledge to Refine Drift Simulations of Early NSS Herring Larvae
Over more than a century, a large quantity of information has been gathered on the life cycle and dynamics of NSS herring (Nakken, 2008). A major part is the identification of several traits during the early larval stages. The main spawning occurs in patches along the continental shelf off the Norwegian coast (Dragesund et al., 1997). The survey data used here for the drift simulations are based on detailed sampling during the main drift phase of early NSS herring larvae along the Norwegian coast. After spawning at the sea floor in February-March , NSS herring larvae hatch about 2-3 weeks later ascending to shallower depths (Blaxter and Ehrlich, 1974) and undertaking diurnal vertical migration (night shallow, day deep) (Munk et al., 1989;Ferreira et al., 2012). The pattern of vertical migration is additionally used to refine drift simulations. During dispersal, larvae stay within water masses of the NCC, where eddies can slow the northward larval displacement often retaining them for a time on banks (Saetre et al., 2002). Hence larval drift was simulated in current fields that should represent the NCC and eddy movements.

Survey Data as Input for the Drift Simulation
Along the continental shelf at <250 m bottom depths NSS herring larvae were collected during and after the main spawning season in March-April . The decision to undertake the north-south extensions was based on information from the spatial extent of the fishery on the spawning grounds, which was modified during the survey when zero larvae were reached. However, as the survey direction was either north to south or vice versa  the boundaries of the larval distribution are not always clearly defined. The transects were perpendicular along latitudes about 28-37 km apart with stations approximately 9.4 km apart. Larvae were collected using two different gears. During day-time double oblique hauls down to 75 m with a GULF III (Gehringer, 1952) (mesh size of 375 µm, mean speed of five knots, nose-cone diameter of 20 cm) were conducted. During night-time, a vertical T-80 ring-net (mouth 80 cm diameter and net 375 µm mesh) down to 150 m replaced the GULF III. The T-80 collected larvae in better condition with similar catch rates compared to the GULF III at night but 10 times lower catch-rates during the day (unpublished data). Station data, important for the use in the simulation, included the location of the station, the time of sampling, the abundance per length class (per 1 mm) collected and the age based on the total length (see below).

Parameterisation of the Drift Simulation
Parameterized station data were used to simulate larval drift within the Lagrangian drift model tool Ichthyop (Lett et al., 2008). This tool is designed for ichthyoplankton dispersal accounting for larval traits like vertical migration behavior, growth, and mortality. For each length class the age was estimated according to a constant temperature influence at 7 • C [mean sea surface temperature at the Møre spawning grounds (Saetre, 2007)]. An approximate growth rate of 0.2 mm per day was implemented (Doyle, 1977). The length of a first-stage yolksac larva was set at 7 mm total length, based on the smallest larvae collected in the surveys. As a test example, we simulated 15 days of backward drift for a larva of 10 mm total length. Each particle was drifted forward until reaching an age of 60 days which approximately mirrors the last preflexion stage 2c from Doyle's (1977) staging system. We restricted the days of simulation according to the youngest larval stages. The reasoning behind this is that older larvae (particularly flexion, postflexion, and later stages) develop an advanced swimming behavior which would increase model result uncertainties (Clark et al., 2005), though assuming the most vulnerable stages being within the very early life stages and thereby probably having the strongest effect on recruitment success (Hjort, 1914).
The starting point of the simulations were fitted with the position of the station where the larvae were collected. The information when a larva was collected was used to account for the most accurate timepoint using particles starting either at night position at 5 m below surface or 40 m below surface when being released during the day (Ferreira et al., 2012). We set the sunrise at 6:00 am and sunset at 9:00 pm -both settings in the spatial-temporal realm being compromises. The starting time was always rounded to a full hour and the end position of a GULF III haul was used as input, whereas the T-80 position reflects a precise position because of the vertical tows.
The hydrodynamic model GLORYS12V1 was used in Ichthyop as hydrodynamic dataset. It is an ocean-eddy resolving reanalysis product with data starting in 1993 freely accessible on the Copernicus Marine Service webpage 1 . The dataset contains a 3D daily mean field on a regular grid approximately 8 × 8 km. The dataset compiled (from global-reanalysis-phy-001-030-daily) comprises daily current fields with 50 vertical sigma levels. Ichthyop has a built-in function reading in the hydrodynamic model component based on the NEMO (Nucleus for European Modeling of the Ocean) 2 framework. The choice of the hydrodynamic dataset was due to usability within Ichthyop and feasibility of simulation duration for a time series covering 24 years of survey observation compromising the relative coarse grid. The advection process was calculated using a 4 th Runge-Kutta numerical scheme. A horizontal dispersion process is implemented in the simulations as the default in Ichthyop (Peliz et al., 2007). The time step used to simulate larval drift was set to 60 s meaning that every minute a new geographical position was estimated. This routine was implemented for each larval size group per station for the back-and forward simulation. New positioning was saved in a NetCDF file for further analyses.

How to Quantify Retention or Drift
Every 4 h the position of a particle from the drift simulations was extracted from the NetCDF files leading to 6 positions per day. These positions were used to calculate the distance (d) travelled per 4 h according to Prasetya et al. (2020): where ϕ is latitude, λ is longitude, R is the earth's radius (mean radius = 6378137 m) and the drift direction (bearing = θ) according to with v larva in m * s −1 (4 h = 14400 s).
To illustrate yearly dispersal patterns 360 • drift plots were produced to identify and summarize the direction and the speed the larvae were drifting. Both the direction and the speed were weighted based on larval abundances observed from the field sampling. Mortality was not included in the model as no specific mortality rates in different areas are known for the investigated periods. Due to limitation of computational power, larval abundances per length class and per station were divided by 10.
Using v and u as Indices for the Velocity and Directional Drift of Fish Larvae A particle (larva) that drifts in water has a direction and a speed. From both direction and speed indices v and u were calculated; v and u are velocities that represent the directional force in the meridional and the zonal direction, respectively. In mathematical terms, v is a velocity only pointing toward a direction on the y-axis while u only on the x-axis. Hence, both can be used as index for either meridional or zonal drift. To calculate v and u a mathematical convention is needed including a conversion of larval drift direction (ldd) to math direction (mathd) according to with θ as mathd converted from degrees to radians. The geosphere package (Hijmans, 2019) was used to calculate the distance and the bearing (direction) per drift segment with the R version 4.0.2 (R Core Team, 2020). A positive v denotes a larva drifting from the south (northward drift), and a negative v a larva from the north (southward drift). A positive u denotes a larva drifting from the west (coastward drift) while a negative u from the east (offshore drift). A yearly index is produced by calculating the mean of v and u for each year and for both the backward v backward and u backward and forward v forward and u forward simulations. The further the mean values of v and u are from zero, the stronger the meridional and zonal drift for the total drift per year.

A Bayesian Approach to Relate Recruitment Success With Drift Indices
The International Council for the Exploration of the Sea (ICES) provides stock assessment indices for NSS herring each year, those of relevance here are recruitment at age 2 years (R age2 ) and spawning-stock biomass (SSB). We used the ratio between R age2 and SSB t−2 (the recruits originating from the SSB 2 years ago)  defined as recruitment success (R age2 /SSB t−2 = R success ), using the assessment undertaken in 2019 which encompasses the 1993 to 2016 year classes (ICES, 2019). The response variable with the larval drift indices were taken as co-variables (v backward + u backward and v forward + u forward ) within a Bayesian inference (Rue et al., 2009). For this the R-INLA (integrated nested Laplace approximation) package within R version 4.0.2 was used.
Data exploration was undertaken according to Zuur et al. (2010). R success index and covariables were a-priori tested for outliers using dot plots indicating no serious issues. Recruitment success in the year 2002 was high compared to other values of the time series. Co-linearity between explanatory variables were tested using the Variance Inflation Factor (VIF < 3). v backward and u backward were more colinear with a VIF < 3 then v forward and u forward with a VIF < 2. Co-variables were standardized (scaled to zero mean and unit variance) and used in two models, to test the drift indices on R success per year, i.e., R success i : where R success_i is µ i with a gamma distribution and a parameter r for the variance, the expected value E(R success_i ) is µ i and µ i an inverse link function (log-link), β 0 the intercept and x i the covariates. A gamma distribution is used as R success is continuous, strictly positive, and right-skewed (positive skewness). Due to the number of 24 years that represent 24 data points in the response variable the two most simple models were constructed with x i representing only drift indices calculated either from the forward or the backward simulations. Two models (model backward , model forward ) were constructed with model backward testing R success i ∼ v backward.std i + u backward.std i and model forward testing R success i ∼ v forward.std i + u forward.std i . The deviance information criterion (DIC) and the Watanabe-Akaike information criterion (WAIC) were used for model selection resulting in the final model which is the full model forward . The final model was validated based on Zuur et al. (2010). Homogeneity was tested visually using dotplots of residuals versus fitted values. Normality was assessed plotting histograms of the residuals. Temporal independence was tested by the autocorrelation function (ACF). Multi-panel scatterplots were used to visually identify relationships between R success and covariables to investigate possible non-linear relationships. Data validation plots can be found in the Supplementary Material (Supplementary Figure 7) and indicate no essential statistical issues.
The final model was interpreted in two ways. First, single effects of the covariables were predicted sketching the model fit for R success based on regular sequences of each range of v forward and u forward . The number of the sequence were based on the number of years from the time series (n = 24). For the final model each sequence either from back-standardized v forward or u forward was combined with corresponding mean backstandardized v forward or u forward values to predict the single effect of the covariables from which the regular sequences were computed. In a second step, the observed versus fitted values of R success were plotted to evaluate overall model performance.

Current Fields Used for Forward and Backward Simulations
Mean current fields from the period 1993-2016 indicate typical circulation patterns along the Norwegian continental shelf area. In the south, between 59 • N and 62 • N, an inner northward flow occurs, branching at around 63 • N north off Møre into two possible drift routes for NSS herring larvae (Figure 1). One branch is closer to the coastline between the coast and the major bank areas with a second flowing northward along the shelf break (∼250 m isobaths). Both branches merge at around 67 • N to the southwest of Lofoten. Up to this point the flow speed is between 0.1 and 0.4 m s −1 . After the merger of both branches at 67 • the flow intensifies frequently exceeding 0.5 m s −1 . At around 70 • N there is a second split into two branches of the NCC, one following the coastline toward the Barents Sea, the other following the shelf break northwards. Stationary mesoscale eddies are formed between the two branches of the NCC, over the major banks (Frøya Bank, Halten Bank, and Traena Bank). Additional non-stationary eddies form off Lofoten and further north where the NCC reaches maximum speeds, which progress offshore. For single years and to visually inspect eddy movement, the reader is referred to video animations ("Video1 -Video6.mp4" found in the Supplementary Material) that illustrates current fields used for larval drift simulations.

Outcome From Forward and Backward Drift Simulations
The larval drift from the backward simulation indicated the highest larval abundances north off Møre as well as on the Halten Bank area which are key spawning grounds of NSS herring (Figure 2). Spawning areas are additionally spread along the continental shelf. While most of the simulated spawning locations were at isobaths <250 m some final positions were also at depths >250 m either in the south of Møre or northwest of Lofoten. At Møre the larval patches drifted mainly northeastwards (Figures 2a versus 2b). From this location onwards many of the larval patches -slightly north off Mørewere transported mainly in the two branches of the NCC taking either the route of the coastal branch or some also along the second branch that follows the continental shelf break. At the widened continental shelf, enhanced larval retention could be observed on the Halten Bank and the Traena Bank (Figure 2b). At the point where the two branches join off Lofoten, larvae were transported farther northwards. Some of the currents and non-stationary mesoscale eddies that form off Lofoten carried larvae offshore. Reaching the shelf break some of the still very young larvae have already drifted into their main nurseries in the Barents Sea, whereas a few larvae followed the second offshore branch further north. During dispersal larvae reached average drift speeds of between 0.1 and 0.5 m s −1 reflecting the typical speed of the currents along the Norwegian coast (Figure 3 and Supplementary Figures 5, 6 360 • drift plots). In all years, the northward drift component was frequent and dominant except for the 3 years 1995, 1997, and 2012, where also a frequent southward drift component was computed (Figures 3a-f)

Larval Drift Indices in Relation to Recruitment Success
The years with a strong southward component that represented enhanced larval retention resulted in years with only low recruitment success (Figure 4). In contrast, the years 2002-2004 with this dominant northward drift component were very successful recruitment years. This is represented by the larval drift indices for the drift along the meridional axis (v backward and v forward ). The only year where an overall negative drift index was computed was from the backward drift simulation for the year 1995. This is the only year where there was a stronger southward than northward larval drift. Hence, in 98% of all years from both backward and forward simulations the northward drift component dominated indicated by positive v backward and v forward drift indices. The meridional drift of larvae was more frequent and stronger than the zonal drift (u forward , u backward ) for both back-and forward simulations (Figures 3, 4). Both mean u drift indices were only positive which means that the overall larval drift per year was inshore. However, years with higher inshore drift coincided with low recruitment years (e.g., 2008, 2010, 2011, 2014) and values closer to zero coincided with years of heightened recruitment (e.g., 1998, 1999, 2016).
From a more detailed view, the success of recruitment was only depicted for some years while the rest of the years could not visually be related to the 360 • drift plots (Supplementary  Figures 5, 6 versus Figure 4a). Although the v forward value in 1996 was small v backward was large and resulted in a higher R success compared to the years 1995, 1997, and 2012 (Figure 4). The  2004, years of low recruitment success coincided with larvae experiencing a strong and frequent northward drift.
The importance of larvae being transported in the meridional and zonal dimension was tested by the full gamma model (model forward ) on recruitment success which was the best model assigned by DIC and WAIC ( Table 1). The fitted and final model was: The meridional drift component v forward indicated an important and positive relationship with R success ( Table 2). The 95% credible intervals of the posterior mean exclude zero which implied that the term is an important factor behind recruitment success.
Simulating v forward from the model forward resembled a continuous increase of R success with stronger v forward (Figure 5a). An increased northward drift until v forward = 0.05 ms −1 apparently constrained recruitment with little enhanced recruitment success and low 95% credible intervals. Higher values allowed for elevated R success by also increasing the 95% credible intervals. That indicated that the potential for higher recruitment is only possible in years where most of the larvae are quickly drifting northwards. Best model predicted u forward values indicate a negative relationship between u forward and R success (Figure 5b). The posterior mean excludes zero values meaning that this index is an important driver of recruitment too, though, in an opposite way; the more u forward is positive the lower is the R success . In other words, when more larvae are transported eastwards toward the coast, the chances of good recruitment falls. So, model forward with only two drift indices as explanatory variables provided an acceptable estimate for recruitment success (Figure 5c). However, the two years 2002 and 2015 were outliers and indicate that larval drift is only one of several factors influencing the highly complex process of NSS herring recruitment.

DISCUSSION
Successful recruitment underlies several episodic and subtle mortality events in early life stages of fish (Houde, 1989). These events can involve very high mortality rates where even small perturbations in environmental conditions can have tremendous effects on the survival of the resulting brood (Houde, 1987).   Larval dispersal is important as it can involve advection to unfavorable environments where larvae may encounter high predation risks and/or low food supply (Hjort, 1926). On the other hand, drift along routes with optimal environmental conditions may enhance larval survival. Knowing about the larval drift patterns that may result in various years of recruitment success is hence an important endeavor to better understand population dynamics and thereby enhance our abilities to manage harvested stocks. Consequently, in short term projections the strength of the Norwegian Costal Current (NCC) should help defining periods with lowered or raised recruitment potential (Skagseth et al., 2015). Such projections can be refined with dedicated survey observations and the generation of larval drift indices. Increasing larval survival rates is dependent on ambient environmental conditions. Here, the adults determine where and when they spawn and thereby ultimately define the conditions their offspring will grow up in Parrish et al. (1981) and Slotte and Fiksen (2000). The spatial limitation, however, is that herring spawning is limited to suitable substrata with appropriate bottom currents (Runnström, 1941;Maravelias et al., 2000). In most of the years, recruitment is low for many herring stocks in the Atlantic and a few exceptional years allow stocks to increase significantly (Trochta et al., 2020). A few weeks in early spring after mature NSS herring have released eggs to the sea floor, larvae will inevitably occur under conditions in which they encounter variable chances of survival. As main transport route, the flow regime regulated by the NCC allows larvae to drift either northwards (Vikebø et al., 2010;Stenevik et al., 2012) or encounter stationary eddies which retain larvae on their way toward the Barents Sea (Saetre et al., 2002). When the NCC is strong, larvae can be carried quickly from spawning grounds northwards (Vikebø et al., 2010). An examination of NSS herring year-class strength and the physical conditions in the NCC (e.g., temperature anomalies, current speed and prevailing wind conditions) revealed that a strong NCC coincides with years in which NSS herring had exceptionally high recruitment success (Skagseth et al., 2015). The findings support the positive relationship between the larval northward advection (strength of the NCC) and elevated recruitment, lending support to Hjort's conclusion that a larval relocation of hundreds of kilometers forms the basis for recruitment success in high-latitudinal waters (Hjort, 1926). This view is also in line with Dragesund's field observations that rapid transport of NSS herring larvae northwards promotes recruitment (Dragesund, 1970). Therefore, it seems that the frequency and the total mean surplus northward drift of larvae control the level of recruitment. This statement is in contradiction of the theory that larval retention in the first months of life may generally be crucial for herring population replenishment (Sinclair and Power, 2015). Especially with the output of the drift simulations in which only one year accounted for stronger retention then northward relocation it seems that the northward advection dominates for early NSS herring larvae. Additionally, the three years in which a stronger southward drift component was observed signify only weak recruitment years.
The NCC is obviously the major transport path for NSS herring larvae to transit from their spawning grounds along the west coast of Norway to their main nursery area in the Barents Sea, but why would the speed of transport be of importance? As temperature rises in early spring the Norwegian Sea becomes more and more accessible for piscivorous feeding during the spring bloom. Juvenile saithe (Pollachius virens), a potential predator of NSS herring larvae, starts schooling at the Norwegian coast in late spring (Husebø et al., 2009). If this schooling coincides with the presence of NSS herring larvae it may lead to an elevated predation on the larvae. To compensate for any possible top-down control by saithe or any other species like Northeast Atlantic mackerel (Scomber scombrus) (Skaret et al., 2015) or Atlantic salmon (Salmo salar) (Utne et al., 2020) it would therefore be beneficial for larvae to escape high predation risks through a rapid transport away from predator aggregations. In this regard, a combination of early spawning by NSS herring and rapid northward transport could combine to enhance larval survival. The earlier drift would enhance the chances of escaping predation but could be disadvantageous due to a lower potential food supply in early spring (Husebø et al., 2009;Slotte et al., 2019). However, the drift speed only seems important in the years before 2004. Most of the years after 2004, coincide with predominant strong northward drift components that failed to produce strong recruitment. Since 2004 regimes seem to have changed in the study area where larvae are susceptible to a fast temperature shift during early development . These altered environmental conditions may also have affected the timing of plankton blooming along the drift route, which normally is delayed toward the north along the Norwegian coast, and resulted in a mis-match with available prey for the larvae . In fact, general shifts in seasonality of Calanus species have been demonstrated over the study period in the Norwegian Sea area, in both coastal and oceanic waters (Strand et al., 2020). Another explanation is that there has been a more general decline of zooplankton production along the entire drift route of the NSS herring larvae as suggested by Toresen et al. (2019). Support for this hypothesis is found in a study looking at the production of Calanus species along the Svinøy section crossing the most important NSS herring spawning grounds, suggesting a 50% decrease for the most oceanic stations, and as much as an 81% decrease for the stations closer to the coast over the period 2000-2011 (Dupont et al., 2017). On the other hand, one cannot exclude the potential of increased predation especially by Northeast Atlantic mackerel that has had a significant northward expansion (Nøttestad et al., 2016), which may overlap and prey directly on herring larvae in the late drift period (Skaret et al., 2015;Allen et al., 2021).
By investigating not only the role of the zonal transport of larvae the model outputs indicate that when less larvae are advected inshore it coincides with heightened recruitment. Nursery grounds in fjords along the Norwegian coast are assumed minor compared to the main nursery grounds in the Barents Sea (Dragesund, 1970;Holst and Slotte, 1998). A drift along and not toward the coast would therefore allow a further transport within the main routes of the NCC. If larvae, on the other hand, would be transported coastwards it may increase chances of larval retention at the coastline or inside the fjords. The drift simulations in this regard could indicate that only a few larvae that were actually collected during the surveys would have been transported beyond the continental shelf that may be low in productivity and may negatively affect larval survival (Norcross and Shaw, 1984). However, the offshore larval transport is area specific. For instance, in areas where productive water masses are mainly observed along regions affected by coastal upwelling provide necessary feeding grounds for offspring (Allain et al., 2007;Tiedemann and Brehmer, 2017). It is hence obvious that in these areas larval retention is essential for survival. Such regions can be observed elsewhere in which fish larvae are suggested to make use of the flow regime to retain themselves onshore (Deschepper et al., 2020;Xing et al., 2020). In the Norwegian Sea, however, the flow regime of the NCC prohibits a strong transport off the continental shelf. Only some currents and eddies that form to the west of Lofoten seem to allow a larval offshore transport well away from shore, but whether this mechanism affects NSS herring recruitment to a significant extent still needs to be investigated.
A valuable add-in to the drift simulations was the actual position, the time, the total length as proxy for age and the larval abundance information from the surveys. With these data, it was possible to produce indices based on real observed spatialtemporal data and indicated years of enhanced larval retention or northward advection. However, even though combining the field observations that provide precision in spatial-temporal dimensions the simulations are marked simplifications of the real environment. The strongest bias may come from the resolution of the hydrographical model used here by probably underestimating small scale processes such as meanders and smaller scale eddies (North et al., 2009;Soufflet et al., 2016). The coarse grid was a choice that could cover the time series examined, the suitability within Ichthyop, accessibility and computational cost. However, drift on scales less than twofold of resolution is cropped that leads to particle trajectories more regular than it would be in real life (Vikebø et al., 2011). This setting has certainly increased the underestimation of the particle spread and probably the overestimation of the branches of the NCC. These limitations would also explain the end positions of the backward simulations that should roughly represent the origin. Besides that most of the origins covered typical spawning areas especially off Møre and the Halten Bank a small fraction of the simulated positions were unrealistic, like the ones at depths >250 m. Thereby, the simulation provides a rough estimate of the larval vertical migration behavior. Further to this, the constant assumed growth rate maximizes uncertainties of the realized drift period until the larvae were sampled in the field. For instance, we used a growth rate of 0.2 mm based on Doyle (1977). However, in other studies a daily growth rate of 0.3-0.6 mm was observed for different larval sizes (Stenevik et al., 1996), whereas in another drift simulation a 0.5 mm daily growth rate was adopted (Vikebø et al., 2011). Hence, the destination as well as origin that we considered in this study should only be noted approximations. Another source of bias may be based on running the drift simulations without region specific mortality rates. In different areas of the study area, larval mortality is most likely variable that would change the frequency of observed drift. For a more complex model it would be therefore necessary to implement spatially resolved mortality rates probably based on predator presence and food availability. A new pilot study was accomplished to use a molecular approach to define predation pressure in the field on NSS herring larvae (Allen et al., 2021). Such work may allow to estimate region specific mortality rates to refine the larval drift model outcomes.
The complexity of fish recruitment variability makes it difficult to assemble all factors that determine recruitment (Houde, 1987). Beside the larval drift there are many other factors that may influence NSS herring recruitment (Garcia et al., 2020). It was beyond the scope of this investigation to test all possible controlling factors; we rather focused on an in-depth investigation into one factor that is a major driver. Investigating recruitment success requires a long time-series that would allow for complex models. Here, the 24 years (24 data points for recruitment success) limits the number of co-variables that can be used in a model, which led us to invoke the "One in ten rule" for linear regression type analyses (e.g., number of observations >10 times the degrees of freedom of the co-variables). However, by developing larval drift indices fitted on survey observations it was possible to link the rapid northward advection of larvae to heightened recruitment years between 1993 and 2016. A strong southward drift component is linked to weak recruitment years. The larval drift indices worked especially well with the years before 2005, since then the changes in the environment (e.g., zooplankton availability and temperature) seemed to have mainly driven the recruitment success (Toresen et al., 2019;Tiedemann et al., 2020). The outcome of this work is intercessional toward larval drift studies in which survey data should be used to refine and calibrate drift simulations. This approach is possible for all planktonic taxa enabling a more precise view on the real drift in planktonic organisms.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
MT, RN, OK, and AS conceived the study. RN and ES collected the data. MT processed and analyzed the data (mapping, drift simulation, and modeling). OK and AS contributed with project resources. MT wrote the manuscript with contributions from RN, OK, AS, and ES. All authors contributed to the article and approved the submitted version.

FUNDING
The work is an outcome of the project RECNOR (Recruitment Dynamics of Commercially Important Fish Species in Changing NE Atlantic Ecosystems) with the IMR project number 14861 funded by the FFA (The Norwegian Fisheries Research Sales Tax System) and additionally funded by the project "Sustainable multi-species harvest from the Norwegian Sea and adjacent ecosystems, " Research Council of Norway with the project number 299554.

ACKNOWLEDGMENTS
The crew members of the various IMR surveys in which the data were collected are gratefully thanked. For data access and curation, we thank Kjell Bakkeplass (IMR-NMD). Alain Zuur and Elena N. Ieno are thanked for an inspiring statistic course in R-INLA that supported the decision for the model choice and presentation.