Exploring drift simulations from ocean circulation experiments: application to cod eggs and larval drift

Drift models are commonly used to study the transport of early life stages of fish and other marine organisms. Various approaches may be applied to examine the distribution and variability of ocean trajectory pathways. In the present study, we compare results using passive Eulerian tracers and Lagrangian float trajectories that are embedded in numerical models. We supplement this analysis by applying an offline model for drift computations. The contrasts in the results from the various configurations are mainly due to differences in drift depth. Simulations were performed using horizontal resolutions of 4 and 0.8 km. The higher-resolution experiment gives somewhat more realistic results for the drift time from Lofoten to the Tromsøflaket bank at the southwestern entrance of the Barents Sea. Furthermore, differences in results between simulation years are much larger than the differences that arise from the choice of model configuration. Climate variability at high latitudes on a multi-decadal time scale is dominated by large interannual variability superimposed on an underlying moderate warming trend. We conclude that a properly configured offline drift model using hourly or 2-hourly results from a simulation with a horizontal resolution of 1 km or finer is the best approach for investigations of trajectory pathways. The flexibility of an offline drift model is also highly advantageous in biological contexts, as it easily allows for a variety of ways in which behavioural characteristics can be parameterized, including descriptions that are defined after the ocean circulation simulation has been executed.


INTRODUCTION
Drift of objects and substances in the ocean may be examined from results of ocean general circulation model (OGCM) simulations. Several approaches may be applied to this end. Objects may be tracked online in simulation time by built-in modules in the circulation model. For substances which are advected and diffused, the same modules that are used for e.g. updating salinity can be applied to general Eulerian tracer fields. Alternatively, particles (Lagrangian floats) may be introduced in an offline model with which positions are updated using stored OGCM output, notably results for ocean currents. Offline trajectory models have been developed for a number of applications, including drift and weathering of oil spills, search and rescue, and transport of ichthyoplankton  and Table 1 therein). Some applications will benefit from additional input such as results for waves and atmospheric conditions; for example, Röhrs et al. (2014) examined the effect of wave-induced drift in a sub-region of the present study's region of interest.
In a recent study, Wagner et al. (2019) examined the representativeness of offline trajectory simulations with respect to results for a passive Eulerian tracer. They concluded that diagnostics such as mean pathways and horizontal distribution were nearly identical for the 2 approaches, while there were some differences in the vertical distributions.
A major advantage over online computations is that offline models are much more efficient, and provide a high degree of flexibility, e.g. for fine-tuning the configuration of trajectory simulations, that cannot be matched by the computationally expensive OGCMs.
In recent years, drift models have become an increasingly popular tool for studying the transport of early life stages (ELSs; eggs and larvae) of fish and other marine organisms (Vikebø et al. 2007, Peck & Hufnagl 2012, Staaterman & Paris 2014. Thus, trajectory modelling has become a tool for advancing knowledge on growth dynamics and distribution of pelagic juvenile offspring. To these ends, both online trajectory models (e.g. Allain et al. 2003Allain et al. , 2007 and offline models (e.g. Strand et al. 2017) have been employed, recording the simulated drift history for several annual spawning events from various spawning grounds. Subsequent examination of the fate of ELSs during drift in conjunction with observations of egg and larval surveys can provide new insights of the space−time variability of fish stock recruitment.
Due to the high abundance and concentration of fish eggs at the spawning grounds, examining drifting objects and advection and diffusion of a substance are both viable approaches. The ocean circulation model that was employed in the present study, the Regional Ocean Modeling System (ROMS), includes modules for online particle tracking as well as an option to add any number of Eulerian tracers to be advected and diffused during run time (Shchepetkin & McWilliams 2005). Both of these approaches are considered in the present study. We are unaware of other studies that have systematically investigated the online tracing of advected and diffused substances as an approach to modelling drift of fish eggs.
The NEA stock of Atlantic cod has increased markedly over the past decades and is currently the largest cod stock in the world (Kjesbu et al. 2014). It sustains significant economically valuable fisheries in several countries, including Norway and Russia (Gullestad et al. 2020). The NEA cod stock is consid-ered to be highly migratory (Bergstad et al. 1987, Robichaud & Rose 2004, with migration distances of 500−2000 km (Opdal & Jørgensen 2015). The population dynamics have typically shown large interannual fluctuations as well as periods of both upward and downward trends (Hylen 2002).
Previous studies have used commercial catch data to investigate variation in spawning ground distribution of NEA cod (e.g. Sundby & Nakken 2008). The spawning of NEA cod takes place in March and April, with a peak around 1 April, but latitudinal differences exist (delayed spawning to the north), and may be related to the seasonal light cycle (Sundby & Bratland 1987, Sundby & Nakken 2008. NEA cod spawn in a number of regions along the coast of Norway (e.g. Sundby & Nakken 2008). Most of the spawning grounds are situated in the general pathway of the Norwegian Coastal Current, and the ELS drift is determined by the strength and dynamics of this current. One major spawning ground is found inside of Vestfjorden, which is separated from the Norwegian Sea by Lofoten, a complex archipelago with islands and narrow straits. In their examination of drift from this spawning ground, Vikebø et al. (2005) found that the dispersion and ambient temperature of drifting larvae were sensitive to the drift depth and only to a lesser degree related to the time of spawning. Recently, Strand et al. (2019) examined the vertical distribution of fish eggs in the same region. They found that the vertical egg distribution determined by wind-induced mixing is modified by vertical velocity shear in the presence of strong horizontal egg concentration gradients. Thus, resolving upper ocean small-scale dynamics is necessary for describing vertical egg distribution, particularly in the early stages.
In their examination, Furevik & Nilsen (2005) found an underlying warming trend in the region of interest for our study. Differences over a period of 30 yr between 15 yr averages from 1958 to 1972 and 1988 to 2002 revealed a rise in sea surface temperature by about 0.2 K or less in the Norwegian Sea. The corresponding warming in most of the Barents Sea was in the range of 0.2−0.5 K (see e.g. their Fig. 19). However, superimposed on the underlying warming trend in the region there is interannual variability of larger amplitudes (1−1.5 K). This variability has been attributed to oscillations in the atmospheric flow, which varies on a time scale of 7−8 yr (Furevik 2001).
This study aimed to explore alternative representations of drift, and to examine how drift patterns and drift environments change in response to interannual variability in the presence of climate change. Further-more, in the framework of the SUSTAIN project, we discuss the application of our results in the context of changing conditions for early life stages of the Northeast Arctic (NEA) stock of Atlantic cod Gadus morhua. SUSTAIN aims to understand how ecosystems can be harvested sustainably by combining the efforts of researchers and end users who are harvesting and managing the ecosystems.
In Section 2, the configuration of the numerical model experiment is described, the various representations of drift are introduced, and the framework for the subsequent analysis is outlined. In Section 3, results on the spatio-temporal variability of drift are given, and the contrasts between the various representations are detailed. In the same section, we investigate how results for drift from Vestfjorden change when we apply a high-resolution configuration that resolves small-scale circulation features and straits in the archipelago. In Section 4 we consider the relevance of results forLagrangian float trajectories and Eu lerian substance advection for drift of fish eggs and larvae. Finally, in Section 5 we discuss strengths and weaknesses of the various methods applied for describing pelagic drift, and make a recommendation for future work on this topic.

Description of the circulation model
ROMS is used to simulate ocean currents, hydrography, passive Eulerian tracers, turbulent mixing parameters and trajectories of Lagrangian floats. ROMS solves the Reynolds-averaged Navier-Stokes equations using terrainfollowing vertical coordinates (Shchepetkin & McWilliams 2005, Haidvogel et al. 2008. The domain and the bottom topography of the model configuration with a 4 km horizontal resolution are displayed in Fig. 1a. The domain covers the Norwegian Sea, the Barents Sea and adjacent ocean regions. The subdivision into analysis regions and their numbering are described in detail in Section S1 in the Supplement at www.int-res.com/articles/suppl/cr01652 _supp.pdf. We also consider a set of larger subdivisions: the Norwegian coast (regions 3−9), the Barents Sea (10−14) and the Norwegian Sea deep basin (region 18−22).
The first mode of the internal Rossby radius of deformation is 3−10 km in the northern Norwegian Sea, with lowest values near the coast (Nurser & Bacon 2014). In the Barents Sea, this deformation scale is 1− 3 km (smallest north of the polar front). Hence, we introduce a complementary configuration with a resolution of 800 m, which covers the ocean off the coast of northern Norway (indicated by a rectangle in Fig. 1a). The bottom topography of this configuration is depicted in Fig. 1b. Hereafter, we will refer to these configurations as 'R4' and 'R0.8', respectively. The R0.8 implementation was defined for the purpose of providing additional information on how smaller-scale circulation features, and the presence of the Lofoten archipelago, impact drift results. The limited size of the R0.8 domain is due to the high computational cost with such a fine horizontal resolution.
In the vertical, the same 36 generalized sigma (scoordinate) levels are defined for both R4 and R0.8 3 Depth (m) The colour coding applied for contouring is displayed by the label bar. Land (dry grid cells) is displayed in gray. The 3 tracer seeding regions are shown in red (M: Møre, Y: Yttersida, V: Vestfjorden). Numbered regions separated by white lines are the analysis regions, detailed in Section S1 in the Supplement (Møre, Yttersida and Vestfjorden are in regions 4, 8 and 7, respectively). The perimeter of the 800 m model domain is drawn as a rectangle in (a). The white dots in (b) inside seeding region V display the initial positions of Lagrangian floats in R0.8. The dashed red line in (b) is the perimeter of the region off Lofoten (to the right) for which results are examined in Section 3.3.3 using the scheme of Song & Haidvogel (1994). Furthermore, stretching of the vertical coordinate that enhances the resolution towards the surface and the bottom was applied, and a minimum simulation depth of 10 m was used in both configurations. The method for computing the horizontal pressure gradient was described by Shchepetkin & McWilliams (2003). For temperature, salinity and Eulerian tracers, the MPDATA advection scheme was applied (Smolarkiewicz 1984). This scheme effectively reduces the numerical dispersion that is associated with many of the second-order advection schemes, thus e.g. preventing negative tracer concentration values. The integration time steps for the barotropic mode were set to 8 and 2 s for the R4 and R0.8 simulations, respectively.
Drift and dispersion are represented in the model integrations applying 2 independent techniques: advective−diffusive Eulerian tracers, and sets of Lagrangian floats that are advected by the local 3dimensional velocity. The results using these techniques were subsequently supplemented by results from offline particle tracking software.

Time-slice experiments
The simulations were performed for a selection of 4 simulation years : 1966, 1980, 1992 and 2006. The simulation periods for R4 and R0.8 were 1 March to 30 September and 1 March to 30 June, respectively. The experiments were initialized on 1 January using results from the SVIM (MET Norway & IMR 2015) archive (Lien et al. 2013), and spun up for a period of 2 mo. As described more in detail in Sections 2.5.1 and 2.4, Lagrangian floats were seeded at selected positions every 2 h during the period 1 March 00:00 h UTC to 30 April 22:00 h UTC, and Eulerian tracer fields were introduced at midnight on 1 March, 1 April and 1 May. The restriction to 4 simulation years was due to limited resources for computations and storage. The selection of years for which simulations were performed was determined in order to represent several decades. Non-processed model results are available from https://thredds.met.no/thredds/ catalog/metusers/arnem-sustainTDS/catalog.html.
Note also that the North Atlantic Oscillation (NAO) index values (Hurrell et al. 2003) during the months leading up to the initialization state (January and February) were mixed (−2.9, −0.7, 1.5 and −0.3 for 1966, 1980, 1992 and 2006, respectively), indicating differences in the winter atmospheric circulation regimes. To examine the strength of the westerlies on a regional scale, we inspected the change in sea level pressure from the Shetlands (1°W, 61°N) to Station M (2°E, 66°N), a position in the central Norwegian Sea. We found that the mean pressure difference in January and February for 1966February for , 1980February for , 1992February for and 2006 were −3.3, 1.6, 8.7 and 4.5 hPa, respectively. Thus, similar variability between the 4 years are seen on larger and regional scales.

Forcing
Daily mean sea level pressure, surface winds, surface air temperatures, surface specific humidity, downward long wave radiation at the surface, downward short wave radiation at the surface and precipitation values were extracted from a regional downscaling of results from the European Reanalysis project (ERA-40, Uppala et al. 2005; ERA interim, Dee et al. 2011). The regional product, NORA10, which has a temporal resolution of 1 h and a spatial resolution of approximately 10 km, is documented by Reistad et al. (2011). Here, we interpolated to the resolutions of the ROMS configurations and updated the atmospheric forcing every 2 h. The atmospheric forcing of the ocean was then implemented by application of the bulk flux algorithm following Fairall et al. (2003).
Daily averages of temperature and salinity profiles, the depth-dependent ocean currents and the surface elevation from the SVIM (2015) archive were applied at the open boundaries of the R4 domain. Also, daily averaged results for sea ice from SVIM were applied on the open boundaries. Details on the configuration for the production of the SVIM archive is given by Lien et al. (2013). Here we note that the ocean conditions at the open boundaries of SVIM were given by results from the Simple Ocean Data Assimilation dataset (Carton & Giese 2008).
In the R4 configuration, freshwater discharge from the coast was represented by monthly climatology values at 19 discharge positions. In R0.8, daily freshwater discharge from 87 rivers were imposed. These discharge data were available from a database constructed by use of the hydrological model HBV (Beldring et al. 2003).
Results from the R4 simulations were applied as lateral boundary conditions in the R0.8 configuration experiments. The atmospheric forcing was derived from the same product as indicated above.

Representations of drift by Eulerian tracers
Tracers were introduced in 3 regions in the model domain, as indicated in red in Fig. 1a. We will refer to these seeding regions as Møre, Yttersida and Vestfjorden (corresponding to regions labelled M, Y, and V in the figure, respectively). The initial vertical distribution of the tracer is shown by the black line in Fig. S2. In the model integration, the tracer was assumed to behave as a dissolved substance with no sources or sinks. Tracer concentration results were stored as daily averages.
The tracer equation solved by the ROMS model is rather complex, and includes options of several turbulence closure schemes. We chose to apply the general length scale model (Umlauf & Burchard 2003) where eddy mixing coefficients are determined prognostically.
With a notable exception, the evolution of the tracer field was computed analogously to salinity. The exception was the treatment of changes in the vertical direction. When diffusive processes act on a passive substance with a given initial distribution, diffusion will act to erode concentration gradients. However, here we wish to implement tracers as a complement to particle representations, and ocean drift tends to be along isopycnals.
Hence, we implemented a simple modification to oppose the tendency towards tracer concentration uniformity in the vertical direction by applying a lift that represents buoyancy. In order to set a value for the lift velocity, we executed a set of test simulations before the production of the model results were performed. Using the same R4 configuration as in the subsequent production simulations, we found that a lift velocity of 0.0013 m s −1 (110 m d −1 ) kept the tendency towards uniformity in check, given the level of mixing in the test simulations.
The principle for implementing a lift velocity as a parameterization of tracer buoyancy is presented in Section S2 in the Supplement. Note that in order to shed light on the impact of this implementation, we inspect a simplified case with constant diffusivity in Section S2. The diffusivity in the model simulations, on the other hand, changes in time and space according to the general length scale model, as stated above.
Tracers were initialized at 3 times in each simulation year (1 March, 1 April, 1 May), in 3 separate regions. Thus, for each year, the evolution of a total of 9 tracer fields was simulated in R4. In R0.8, only the evolution of the tracer fields initialized in Vestfjorden on 1 April are considered.
The tracer results from R0.8 were aggregated onto the coarser R4 grid on 1 June, and then introduced as initial conditions in the R4 simulations for the period 1 June to 30 September. The restriction to a 2 mo period with R0.8 was imposed in order to retain virtually all tracer mass inside of the R0.8 domain at the time of re-initializing tracer concentrations in the R4 simulations. Simulations using this hybrid configuration will hereafter be referred to as the 'combined experiment'.

Representations of drift by Lagrangian floats
Using floats to represent drifting particles has several strengths when it comes to flexibility in examinations of results. The particles can be assigned with a non-uniform distribution of a property. Particles seeded at the same geographical position, but at different times may e.g. be taken to represent a timedependent number of fish eggs. Similarly, particles seeded at the same time may be assigned to represent a spatially varying initial distribution of fish eggs.
For simulation of drift by floats, the spatio-temporal variability at seeding may be considered a posteriori. Nevertheless, in the present study such variability is not considered. Since the environmental history (e.g. ambient temperature) of each particle may be recorded, parameterizations of particle-specific properties (e.g. mortality) may also be taken into account a posteriori.

Online floats
Unlike a tracer, which is a grid cell property in the model (see Section 2.4), online floats represent Lagrangian particles with given longitude−latitude− depth positions that change with time. The Lagrangian floats drift with the 3-dimensional velocity fields as given by the results from the circulation model, including results for the vertical velocity, with no further modifications such as a lift velocity. Floats were released into the model simulation every 2 h during the first 2 mo of the experiments, i.e. during the period 1 March 00:00 h UTC to 30 April 22:00 h UTC.
In the R4 simulations, floats were seeded on 24 positions (9 positions inside each of the seeding regions Møre and Yttersida, and 6 positions in Vestfjorden). At all positions, floats were released at 3 depths (5, 15 and 35 m). Thus, 72 floats were released every 2 h, giving a float count of 52 704 for the 2 mo release period for each year, i.e. the total number of simulated float trajectories for the 4 simulation years was 210 816.
In the R0.8 simulations, floats were released at the same depths at 12 positions inside seeding region Vestfjorden. These positions are shown as white bullets in Fig. 1b. Complementary 'combined experiment' results were obtained by inspecting float results from R0.8 after 2 mo of drifting, and positioning them in the R4 grid at that time, and then continuing the trajectory simulations until 30 September. The restriction to a 2 mo period with R0.8 was imposed in order to retain all floats inside of the R0.8 domain at the time of re-initialization.
Trajectory positions and depths were stored every hour, along with complementary information of the float environment (e.g. temperature and salinity along the trajectory).

Offline floats
As a supplement to the online simulations, we also applied the offline drift software 'OpenDrift' . OpenDrift is a generic framework for trajectory simulations which has been written to facilitate advection of particles chosen from a set with wide-ranging properties, e.g. representing oil, fish eggs and capsized boats. For the present study, we applied a simple choice of buoyancy-modified vertically mixed particles from the class OpenDrift3D as given e.g. in Fig. 6 of Dagestad et al. (2018) for baseline offline simulations. Particles were seeded in the same positions and depths and at the same times as the online floats (see Section 2.5.1 for details). In this baseline, we applied a forward Euler scheme. The internal time step in OpenDrift was set to 30 min. For the implementation of buoyancy effects, the same lift velocity as for tracer concentration simulations was applied (0.0013 m s −1 ); see also the discussion in Section 2.4 above, and the conceptual and simplified presentation in Section S2 in the Supplement.
The OpenDrift trajectory software is highly flexible, as drift may be simulated using a range of available computation modules. In addition to the baseline described above, we performed an extensive set of simulations using alternative modules and config-urations. The baseline was configured so that stranding at the coast or at the bottom was denied by returning stranded particles to their pre-stranded position prior to executing the next integration time step. We performed control simulations where stranded particles became immobile, and report on the results in Section 3.3.3. We also checked the use of the forward Euler scheme by comparing with results obtained when applying a fourth-order Runge-Kutta scheme. With one exception, which we mention in Section 3.3.3, the results became nearly identical to the OpenDrift baseline experiment.
Furthermore, we supplemented the simulations with alternatives where we either modified the lift velocity or applied a module for drift of fish eggs and larvae , both of which are approaches that are more in line with the formulations used in the existing literature on drift of cod eggs and larvae. Results from these simulations are discussed in Section 4 (see Fig. 2).

Subsetting results for online floats
While batches of online floats in R4 are released every 2 h during a 2 mo period for each simulation year (see Section 2.5.1), tracer releases are introduced 3 times each year, at the start, middle and end of the 2 mo period. Thus, when comparing results from R4 for tracers and floats in Section 3 below, we limit the online float releases that we include in the analyses to those released during the initial, middle and final week of the 2 mo period.
From the R0.8 simulations, we have tracer concentration results from the 1 April initialization, so for comparisons of R0.8 results, we limit the analyses to floats that are released in the middle week.
As explained in Section 2.4, a lift is introduced in the computation of tracer concentration, to compensate for the tendency towards a uniform distribution in depth due to diffusion. For this purpose, we applied a vertical velocity of 0.0013 m s −1 (110 m d −1 ) upward. Integration of tracer concentration results from this configuration revealed that on average, only 0.2% of the total tracer mass was displaced to depths >100 m.
The vertical velocity fields in ROMS are determined diagnostically from the continuity equation, which leads to noisy results. The numerical solution technique thus introduces numerical dispersion in the vertical distribution of floats, akin to the effect of diffusion that we discussed in Section 2.4. So ideally a lift velocity should have been applied when updat-ing the vertical position of the online floats, but this was not implemented at run time. When inspecting the results for online floats, we found that about half of the float trajectories reached depths in excess of 100 m at some time.
As a consequence, we discarded all results for online floats which recorded vertical positions beneath the 100 m level. Such a restriction is also in line with the vertical profiles from theory as well as from observations of fish eggs, e.g. Solemdal & Sundby (1981) and Sundby (1997).
The reduction in the full set of online floats after imposing the depth criterion is not uniform in time or space. To compensate for the contrasting distributions, weights are assigned to the floats in the reduced set so that at each seeding time and for each seeding region, the same total weight is applied.

RESULTS
We examined how average lifetime depths in the representations evolve as a function of drift time. The results for each of the seeding regions are displayed in Fig. 2. For floats, this was achieved by averaging the vertical positions of the trajectories from the time of seeding to the age as given by the values on the horizonal axes. Recall that floats are initialized 2hourly over a period of 2 mo, and thus reach a certain age over a period that also spans 2 mo. For tracers, the averages of the daily vertical centre of gravity (sometimes referred to as the centre of mass) were used. Furthermore, these mean depths were computed by averaging results for all simulation years. We found that the lifetime mean depths as functions of age for tracers and for floats reach a stable level after about 2 wk.
We note from Fig. 2c that in R4 the representation by tracers on average is elevated by 5−10 m above the online floats (solid lines in Fig. 2a), owing to the buoyancy parameterization described in Section 2.4. Moreover, even though the same lift velocity was applied for the results displayed in Fig. 2b and c, the offline floats (panel b) have somewhat deeper drift depths. It is worth noting that there are conceptual differences between Eulerian tracers and Lagrangian floats, e.g. the parameterization of mixing effects for the 2 implementations gives rise to differences in the vertical position.
We also found that in the longer perspective, a lift velocity of 0.0013 m s −1 is somewhat on the high side, since this lift raises the centre of gravity by several metres above its initial state. We have supplemented our results with offline trajectory modelling applying   (Solemdal & Sundby 1981, Sundby 1983). The effect of a reduction of the lift velocity from 0.0013 to 0.001 m s −1 is a deepening of the average drift of about 3.5 m (full vs. dashed black lines in Fig. 2b). For R0.8, the drift of online floats takes place at about the same vertical level as in R4. However, the tracer and offline floats, both with buoyancy effects parameterized, are lifted to higher vertical levels in R0.8 than in the corresponding tracer results from R4. These contrasts must be kept in mind when evaluating the simulation results for the various representations of drift.
Furthermore, we have integrated the results for floats and tracers over the analysis target regions and the larger subdivisions that are described in Section S1. The results are discussed in detail below, for the various seeding regions and for the set of model configurations.

Drift from Møre
The largest changes in distribution results in Table 1 between the 4 simulation years were found for the drift from Møre. Here, the average concentration in the Barents Sea regions during 1966 is larger than for the other simulation years by a factor of 3−10, for both tracer results and for floats (not shown). One important aspect in this context is the drift time northward along the Norwegian coast to the Barents Sea. We examined the drift time (age) of tracers and floats at which their fractions first exceed 2% in the various analysis regions, for the full set of tracers and for all float release times. With this definition, we found that the drift times to analysis regions 8 and 9 (Lofoten/Vesterålen and Troms coast north, respectively) were 2−6 wk shorter in 1966 than in the other years. We recall from Section 2.2 that 1966 had a much lower January−February NAO index than the other 3 years. Hence, the contrasting pre-conditioning for the ocean circulation in spring and summer in the simulation years may be of relevance for the rapid drift in 1966.
To further illustrate the contrast between results for 1966 and other simulation years, the temporal evolution of float fractions from Møre in the Barents Sea regions is displayed in Fig. 3. Shown are results based on floats that are seeded during the first week of March, for all 4 years. The rapid drift in 1966 is evident as we note that in simulation year 1966, the Barents Sea float fraction exceeds 2% almost a month earlier than in 1992, and more than 2 mo earlier than in 1980 and 2006.
Inspecting the time of 2% exceedance for tracer substance from Møre to analysis regions 8 and 9, we find that drift times are 1−3 wk shorter than the corresponding drift time for floats. The slightly swifter drift of tracers than floats is also seen in most annual results (but not all) for drift from Yttersida and Vestfjorden. This contrast is due to currents being somewhat stronger for the tracers which on average are elevated by 5−10 m above the floats (Fig. 2a, Results for floats seeded during 3 selected weeks from the 2 mo release period are provided in Table 2. For this set of weeks, the drift time to the 2 northernmost regions among those that have been tabulated (9 and 10) is shortest in 1966. Long drift times are associated with high temperatures, as expected due to the seasonal evolution of ocean temperatures in the upper 100 m. The swiftly drifting floats in 1966 generally have lower average temperatures when their relative concentration reaches 2% in 'Barents Sea SW'. An extreme case is the initial release week in 1966, when cold conditions led to mean temperatures below 5°C for floats in the Barents Sea.

Drift from Yttersida
From Table 1, we note that about half of the drift time in the 5 mo period is spent in the Barents Sea in the R4 simulations. Fig. 4c,d reveals that while the largest concentrations occur in 'Barents Sea SW', a notable amount (about 20%) of the drift from Yttersida takes place inside of the 'Barents Sea central' region.
Given the strong currents in the region from Yttersida to the Barents Sea that exist in both model results and observations (Melsom & Gusdal 2015), it is no surprise that the drift to 'Barents Sea SW' (region 10) is rapid. Within a month, the regional concentrations reach 5% in all simulation years, and within 2 mo, the concentrations reach 20%. The conditions that tracers and floats encounter as they move from Yttersida to the Barents Sea are colder by about 1−2 K in simulation years 1966 and 1980 when compared to 1992 and 2006, except for the final release week in 1980 (Table 3). Also, the temperature conditions are fairly uniform along trajectories, since the temporal averages change by 0.5 K or less. Exceptions to this can be noted for a few cases when the drift time exceeds 2 mo, and for the last release week in 2006.
The number of floats which are discarded due to the depth restriction is by far greater for those seeded in Yttersida than for the other seeding regions. In this context, one should note that Yttersida is closer to the deep basin than the other regions. An inspection of the trajectory results reveals that while the full set of trajectories from Møre and Vestfjorden are positioned over ocean depths > 250 m about 20% of the time, the corresponding value for trajectories that originate in Yttersida is 40%.

R4 results
The seeding regions of Yttersida and Vestfjorden are positioned with a modest degree of spatial separation: the distance between their centres is less than 100 km. Nevertheless, when we compare the regional representation of R4 results for floats (e.g. Fig. 4d,f) the distributions are noticeably different. The contrasts between the corresponding results for tracer distributions (Fig. 4c,e) are even more pronounced. The disparity between results for drift from Yttersida and Vestfjorden is also evident in Table 1.
Furthermore, we examined the time before the horizontal centres of gravity became displaced by more than 200 km from their initial positions.  Table 4. For all years, this displacement occurs at a much later stage for the Vestfjorden results, for both tracers and floats. The contrasting results are likely affected by the positioning of Yttersida close to the shelf break (Fig. 1a), where currents are strong. Furthermore, the seeding regions Yttersida and Vestfjorden are separated by the Lofoten archipelago, with Vestfjorden on the inside of Lofoten. The archipelago coastline is poorly resolved in R4, and its land−sea mask represents a peninsula as all  straits are closed. The results such as the time taken to reach the Barents Sea may well be substantially impacted by the spatial resolution. In order to examine how the simulated drift is affected by the resolution, the R0.8 experiment was introduced.

R4 vs. R0.8
First we inspect results from the combination experiment as defined in the last paragraph of Section 2.4, and compare them to R4 results for drift from Vestfjorden. The time for 200 km separation from the initial positions listed in Table 4 reveals that the results for tracers are only modestly impacted by the configuration choice. However, for the floats, the time interval to reach 200 km separation distance is longer in the combined experiment, and the contrasts to the results for seeding region Yttersida are further inflated. To exemplify, consider the results for Vestfjorden in 1966. The time at which the tracer gravity centre is at a distance of 200 km from the seeding position is 60 and 62 d with the R4 and the combined configuration, respectively, a difference of only 2 d. The corresponding times for 200 km separation for floats in the 2 configurations are 71 and 96 d, respectively, i.e. a difference of 25 d.
The shorter drift time for floats in the R4 configuration is consistent with the results for the drift time when the fraction of floats that reach 'Barents Sea SW' exceeds 5%, although for one year (1992), the drift time is nearly identical (98 and 99 d).
Next, we examine how the lengths of the online float trajectories, as measured along their tracks, are affected by resolution. We supplement the raw results for trajectory lengths by application of a temporal box filter with an integration period set to 2 periods of the M 2 tide component (24.84 h). This simple approach also effectively filters out motion for smaller, but nonnegligible daily tidal constituents. Details and results for a sample trajectory are given in Section S3 in the Supplement.
Statistical distributions for the set of trajectory length results are shown in Fig. 5 after 2 mo drift. The raw trajectories (Fig. 5a) in R0.8 are longer, thus with higher average speeds, than the float trajectories in R4. Furthermore, the same inference is made for the lengths of the detided trajectories (Fig. 5b), so this difference between configurations is not due to differences in the tidal motion. This result may initially seem to be at odds with the swift separation in R4 that is evident from Table 4. However, even though the float speed is higher in R0.8, the trajectories spend a longer time in analysis region 7 (Vestfjorden) in R0.8 (65% of the total drift time) than in R4 (57%). Hence, R0.8, which resolves circulation features with smaller horizontal scales than the R4 configuration, exhibits a more effective retention inside Vestfjorden than the R4 configuration.

R0.8 results
We now consider results from the set of R0.8 experiments, including the offline float simulations with the OpenDrift model (see Section 2.5.2 for details). Note that these simulations are limited to the period 1 March to 30 June for each of the 4 simulation years.
We first consider the mean vertical position of tracer mass and floats for the releases that correspond to the 1 April tracer release. The results are depicted in Fig. 2d. The ROMS floats and tracer mass are centred at depths of about 20 and 4 m, respectively. This is a larger contrast than we found with the R4 configuration experiments, where the corresponding depths were approximately 13 and 6 m, respectively (black lines in Fig. 2a,c). Furthermore, the mean depths from the OpenDrift offline simulations closely match the corresponding results for the tracer mass depth, as expected from the parameterization of buoyancy (see Section 2.5.2).
Next, we examine float trajectories that enter the region offshore of the Lofoten archipelago, given as the region to the right (north) of the dashed line in Fig. 1b. The drift time at which floats enter this region is displayed in Fig. 6. Offline floats (Open-Drift) that passed to the south of the Lofoten archipelago arrived slightly earlier than online floats (ROMS). This is due to the difference in the vertical  Table 4. Number of drift days before the centre of gravity displacement exceeded 200 km from its initial position. These results are based on simulations of Eulerian tracers initialized on 1 April and Lagrangian floats released during a week centred on that day. Columns Y and V give results based on drift from seeding regions Yttersida and Vestfjorden, respectively. The 'Combo expt' columns give results for the combined experiment as defined in Section 2.4 (only results from Vestfjorden are shown for the combined experiment) position (Fig. 2), since currents generally have higher speeds near the surface. Furthermore, the total float fraction off Lofoten after 60 d of drift is nearly the same with online and offline models (45 vs. 42%, respectively), and the corresponding fractions that leave Vestfjorden through the straits are also similar in the 2 cases (4 vs. 6%). Recall from Section 2.5.2 that stranding of Open-Drift particles was suppressed. We also examined a set of complementary simulations in which trajectories were terminated upon stranding. This dramatically changed the results, and the passage of offline floats through straits was nearly eliminated, as only a fraction of 0.05−0.18% (for Euler and Runge-Kutta, respectively) then passed through the straits. Note that theoretically, currents do not cross boundaries (land, bottom) in our configuration, which is why there is no stranding of online floats. In the present implementation for computations of offline float trajectories in OpenDrift, the situation is different; vertical motion is parameterized by buoyancy-modified mixing, and horizontal velocity updates are available with a limited resolution of 2 h.

Synthesis of drift results
Here, we compare the distribution of tracer mass and floats from a variety of configurations, after 60 d drift from the Vestfjorden seeding region. The fractional distributions in the various analysis regions are displayed in Fig. 7. The differences in the results from the various configurations are modest. When examining the results from the R4 experiment (Fig. 7a,c), the downstream concentrations are slightly higher for tracers than for online floats. This is related to a somewhat faster drift of tracer substance when compared with floats, similarly to the corresponding results from seeding region Møre (Section 3.1). The concentration of floats shown in Fig. 7c (R4) and Fig. 7d (R0.8) is indicative of larger retention in the latter, with larger concentrations inside Vestfjorden.
Results for offline floats are slightly different, with 46% of the floats retained inside Vestfjorden after 60 d with the R4 configuration results for OpenDrift simulations. The corresponding fraction with the R0.8 configuration results is slightly lower (42%). Hence, the retention rate may be sensitive to the drift depth, which is deeper for online floats (e.g. Fig. 2d).
Finally, the results from the R0.8 online and offline simulations (Fig. 7d,f, respectively) exhibit a high degree of similarity. However, there is a slight con-  (Fig. 2d), which is associated with slightly lower drift speed.

Interannual variability
Interannual variability is a dominating signal of climate change in water mass properties in the eastern Norwegian Sea and the Barents Sea, for time scales of about 50 yr or less (Furevik 2001, Furevik & Nilsen 2005. By selecting years with contrasting NAO indices, we have found that the interannual atmospheric variability gives rise to differences in drift pathways and the physical environment for fish eggs and larvae. The quality of the drift results depends decisively on the accuracy of the model results for ocean currents. Unfortunately, the amount of direct and indirect observations of ocean currents from the 4 simulation years is much too sparse for a proper evaluation of the quality of model results. In Section S4 in the Supplement, we present results from model validation of salinity and temperature in the upper 100 m. This validation reveals the high quality of model results for temperature along the Norwegian Coast as well as in the Barents Sea, both with respect to spatial and spatio-temporal (interannual) variability.
Samples that display interannual variability in drift pathways from the simulations with the R4 configuration are presented in Section S5 in the Supplement. The drift pathways in our results exhibit a notable variability, e.g. in the distribution along the coast, and the degree to which drift is displaced off the coast. Moreover, there is a distinct variability be tween years in the retention inside Vestfjorden of Lag ran gian floats and Eulerian tracers (Fig. S8a−d). See Section S5 for further examples of interannual variability that the small sample of trajectories illustrates. Bjørke & Sundby (1987) examined survey data for the geographic distribution of the smallest post-larval cod. They found that the fractions of the total post-larval cod on the deep side of the eastern continental slope of the Norwegian Sea (their regions 1 and 2), were 0.03, 0.35 and 0.18 for the years 1983, 1984 and 1985, respectively (their Table 1). In 1988, the fraction of pelagic juveniles in this deep-sea region was estimated at 0.35 by Suthers & Sundby (1993). Using 38 yr of ocean current results from the SVIM archive as input to offline trajectory computations with LADIM, Strand et al. (2017) found that an average fraction of 0.115 of the trajectories from 10 spawning regions was situated in the deep ocean on 16 September (increasing to 0.147 when weights for spawning intensity were introduced). In the present study, we identify target regions 18 and 20 as largely overlapping with the deep-sea regions of Bjørke & Sundby (1987). When we analyse the results from our 4 years of online trajectory simulations for 16 September, the average fraction of floats in this deep-sea region is 0.12. This is close to the result given by Strand et al. (2017) (regions 1, 4 and 5 in their Fig.  2A), although the fraction they reported is slightly higher. As one would expect from e.g. Fig. 4, the corresponding fraction of tracer mass is higher. For the 4 years under consideration, the average tracer mass fraction in this deep-sea region is 0.25 on 16 September.

RELEVANCE FOR DRIFT OF ELS OF FISH
In an earlier study of transport of cod larvae from Lofoten, Ådlandsvik & Sundby (1994) found that the transport in their model was too fast in the initial phase. While surveys indicate that it takes at least 80 d for the major part of the population to reach Tromsøflaket (southwest in our region 10), their simulation results suggest a drift time of about 50 d. They suggested that the rapid drift may be related to unresolved coastal and bottom topography features in their configuration with a resolution of 20 km. The spawning area which they examined is similar to a combination of regions Vestfjorden and Yttersida in our study. Hence, a corresponding quantification of drift time in our results will depend on the relative weights assigned to drift from the 2 regions. If we disregard results for the outlier year 1966, we find that the drift time from Yttersida is about 50 d, as in Ådlandsvik & Sundby (1994). However, the drift time from Vestfjorden in the results from R4 is 56−89 d for tracers and 74−98 d for floats (Table 5), i.e. close to the drift time inferred from survey data.
Furthermore, as noted in Section 3.3, the drift time increases when the horizontal resolution be comes finer. In the case of the combined experiment, the drift time to 'Barents Sea SW' generally increases with an order of 1 wk (Tables 4 & 5). Based on the analysis in Section 3.3, this was due to increased variance in the results with the R0.8 configuration, likely due to resolution of smaller-scale processes. This is in line with the suggestion of Ådlandsvik & Sundby (1994) that adopting a finer resolution in the hydrodynamic model would lead to an improvement of their results. With a horizontal resolution close to the present R4 configuration, Vikebø et al. (2005) found that retention inside Vestfjorden can change from one year to the next.
As stated in Section 2.2, the simulation period with R4 is 1 March to 30 September. The start of this simulation period was set according to early spawning of NEA cod. Further, the 0-group is still fully pelagic in August, but becomes more associated with midwater layers and undertaking diurnal vertical migrations. By October−November, the juveniles have usually become predominantly demersal . Hence, we chose to terminate our simulations at the end of September. Note also that the distribution of spawning changes in response to phase changes of the Atlantic Multidecadal Oscillation (Sundby & Nakken 2008). In the present context, such variability can be taken into account in the processing of model results, as outlined in Section 2.5.
The 3 spawning areas of NEA cod that defined our seeding regions have all been major regions of spawning at some point. Nevertheless, this is not a complete set of NEA cod spawning areas. For a more comprehensive review of spawning areas, readers are referred to e.g. Fig. 1 in Langangen et al. (2019).
As is seen from Fig. 2, the various representations span a range of depths. Thus, the results for e.g. tracer and floats from R4 can be taken to indicate bounds for the drift distribution. Note, however, that even though the drift depths for tracers are closer to the surface than the drift of online floats, the geographical distribution as displayed in Fig. 4 exhibits relatively small differences (panels a,c,e vs. b,d,f, respectively; also a,c vs. b,d in Fig. 7).
On the other hand, as revealed e.g. by Table 5 and Fig. 3, contrasts are large between the 4 simulation years. Thus, our results suggest that there have been substantial changes in the drift conditions for eggs and larvae. Hence, interannual variability and longterm changes in oceanographic conditions are highly relevant when assessing the potential benefits in spawning time and location of NEA cod (Opdal et al. 2011, Langangen et al. 2016. We performed auxillary simulations of offline drift with a fish egg/larvae module in OpenDrift, to shed some light on the impact of changes from the passive transport that have been considered so far, to a representation of ELS of cod. Simulations were performed using results from the R0.8 configuration, and the average lifetime drift depth from Vestfjorden is displayed by the blue line in Fig. 2d. The level of drift with this implementation is deeper than the regular simulation for offline floats (solid black line). This is in line with the discussion of the magnitude of the lift velocity in Section 3, where we presented results using lift velocity based on the buoyancy of cod eggs from Solemdal & Sundby (1981).
Furthermore, we need to stress that in this study we did not consider a complementary aspect for survival success of ELS, namely the availability of food. Skogen et al. (2007) found that primary production in the northeastern Norwegian Sea is negatively correlated with NAO. We also did not explicitly account for spatial variation in natural mortality (Langangen et al. 2014), caused by e.g. starvation or predation.

CONCLUDING REMARKS
We have examined and compared results from various approaches to drift modelling. Results using Eulerian tracers and online Lagrangian floats are similar, and differences can be attributed to effects of the configuration of drift simulations. The contrasting distributions of trajectories from the various approaches are manifested in differences in the drift depth.
The advantage of using Eulerian tracers is that a more complete initial state can be prescribed, leading to results that represent evolution in a continuum rather than with a limited set of individual trajectories. This way, results are more complete in the sense that all resolved circulation features have an impact on the distribution of the substance that is drifting.
On the other hand, online Lagrangian floats can be specified with a degree of flexibility that in practice cannot be matched by Eulerian tracer representation of drift. Results for a particle represent a trajectory from a known starting position. Each time step at  Table 5. Age in days when the fractions of Lagrangian floats and tracer mass seeded in Vestfjorden first exceeded 0.05 in selected target regions. Region specifications are provided in Fig. S1. Tabulated tracer ages are computed by integrating all 3 Eulerian tracer releases, while drift times (ages) for floats are given for all release times. 'combo expt' is the combined experiment (defined in Section 2.4) which a Eulerian tracer is seeded requires simulation of a full 3-dimensional field, which in R4 corresponds to about 10 7 positions. Each Lagrangian float corresponds to only 1 position at each time step. However, note that each float may be taken to represent different initial conditions, as explained in Section 2.5. In the present study, we seeded sets of Lagrangian floats every 2 h during a 2 mo period. Due to large computational demands, Eulerian tracers were seeded only 3 times. Online floats are economical not only when it comes to computing time but also when disk space required to store results is considered. The third approach that we have examined, drift of offline Lagrangian floats, is by far the most flexible strategy for drift simulation: no decision regarding seeding times and positions must be made prior to the execution of the ocean circulation experiment. Offline computation of drift is much less computationally demanding, so e.g. discrepancies in one configuration for drift simulation can rapidly be corrected in a follow-up. The disadvantage with offline trajectory modelling is that this approach is limited by the amount of information that is available from the circulation experiment that provides input to the drift simulation. Updating trajectories is susceptible to inaccuracies from interpolation since the full temporal evolution of the ocean current field will in practice not be available. If not properly addressed, these shortcomings can lead to excessive stranding of trajectories. From their examination of Eulerian tracer concentrations and offline particle drift in a deep ocean setting, Wagner et al. (2019) suggested that numerical dispersion and their use of daily mean currents applied for particle advection are possible reasons for discrepancies between results from the 2 approaches.
A supplementary method to gain insight on trajectory pathways from offline simulations is to perform experiments by adopting backward-time mode. However, in the present case of non-linear motion in a turbulent, diffusive ocean, there are issues with irreversible processes that make backtracking problematic (Batchelder 2006). Hence, we have refrained from adopting backward-time mode offline simulations in the present study, but we acknowledge that such an approach can be useful if results are interpreted with the limitations of its validity in mind.
Aspects of configuration for drift experiments that have not been considered in this study include Stokes drift due to wave motion (Röhrs et al. 2014) and dispersion effects on Lagrangian float trajectories. The latter may be implemented by e.g. a random walk parameterization both in online and offline computations, and again the flexibility of the offline approach is relevant. Barring the use of a much more sophisticated model system in which wave and ocean circulation simulation is run as a coupled system (e.g. Warner et al. 2010), effects of wave motion can only be implemented with an offline configuration. Note, however, that even in an offline model there are challenges when combining results from separate simulations of waves and ocean circulation, since generally all of the momentum transfer from the atmosphere is consumed in ocean circulation experiments. In principle, no momentum is then available for generation of waves.
In Section 3.3.3, we found that retaining information as 2-hourly output from ocean circulation simulations that include tidal motion was sufficient for realistically reproducing the temporal evolution of the distribution of trajectories as they exit the Vestfjorden region in an offline trajectory model. In Section 4, we noted that drift time from Vestfjorden to Tromsøflaket was more realistic in results derived from the higher-resolution R0.8 configuration. Hence, our recommendation from the present study is to perform an eddy-resolving ocean circulation experiment, store results every hour or every 2 h, and apply these results with an offline drift model that has a high degree of configuration flexibility.