Sensitivity of modelled North Sea cod larvae transport to vertical behaviour, ocean model resolution and interannual variation in ocean dynamics

Abstract

Transport with ocean currents affects the spatial distribution and survival of fish eggs and larvae and thereby population connectivity. Biophysical models are commonly used to understand these dynamics. Advancements such as implementing vertical swimming behaviour and higher resolution ocean circulation models are known to improve model performance, however, the relative importance of vertical behaviour vs. ocean model resolution is elusive. Here, we use North Sea cod (Gadus morhua) as a case study to assess how vertical movement, ocean model resolution and interannual variation in ocean dynamics influence drift patterns and population connectivity. We couple a fine (1.6 km, 3 h) and coarser (4 km, 24 h) ocean model to an individual-based model for cod eggs and larvae, and compare simulations with and without vertical movement of eggs and larvae. The results are moderately influenced by vertical movement and ocean model resolution but differ substantially between years. While ocean model resolution is consistently more influential than vertical movement, the effect of vertical movement strongly depends on the spatiotemporal scale of the analyses. This study highlights which aspects of biophysical modelling of connectivity that most critically affect the results, allowing better investing computational resources and proposing goal-based guidelines for future studies.

Atlantic cod, biophysical model, Gadus morhua, interannual variability, larval drift, model sensitivity, North Sea, particle tracking, connectivity, vertical movement

Introduction

For many fish species, such as Atlantic Cod (Gadus morhua), the success of a year class is strongly dependent on survival during the first year (Cushing, 1990), when a long period is spent as planktonic larvae. Quantifying drift patterns is important to understand drivers behind variable larval survival, such as temporal or spatial overlap with food (the “match/mismatch hypothesis”, Hjort, 1914; Cushing, 1990), favourable retention (the “member/vagrant hypothesis”, Sinclair and Iles, 1989), or advection to nursing grounds (the “migration triangle hypothesis”, Harden Jones, 1968). Moreover, drift patterns can play a role in maintaining population structure (e.g. Heath et al., 2008), with important implications for management of harvested populations (Heath et al., 2014).

Biophysical modelling is a common tool to investigate drift patterns of fish eggs and larvae. During the past decades, we have seen improvements in the resolution of ocean circulation models used in biophysical models, and advanced our knowledge about larval behaviour, which is incorporated into increasingly complex individual-based models (Staaterman and Paris, 2014). Numerous studies have shown how implementing vertical movement influences drift patterns and may increase retention (e.g. Cowen et al., 2006; Fiksen et al., 2007; Vikebø et al., 2007; Ospina-Alvarez et al., 2012). Ocean model resolution can also significantly influence drift patterns (Huret et al., 2007; Qin et al., 2014). For example, Putman and He (2013) improved the correspondence between modelled and observed drift trajectories of juvenile sea turtles by using a high-resolution ocean model relative to coarser products. Still, few modelling studies assess results’ sensitivity to ocean model resolution (but see Hufnagl et al., 2016), and, to our knowledge, the relative importance of ocean model resolution vs. implementing vertical movement has not been investigated.

Faced with the computational cost of running high-resolution ocean models and complex individual-based larval modules, researchers typically have to compromise between computational cost and model realism. Similarly, increasing model complexity may constrain the number of years that can be simulated. To help investing computational costs more efficiently, we here compare the relative sensitivity of drift model results to (i) implementing realistic vertical movement of eggs and larvae, (ii) improving ocean model resolution, and (iii) inter-annual variability in ocean dynamics, using Atlantic cod in the North Sea as a case study.

Albeit currently being managed as one stock, North Sea cod can be divided into two genetically distinct units (Figure 1): Viking in the northeast and South (Dogger) in the south and west, likely being separated by a combination of adult migration and larval drift (Heath et al., 2008, 2014; Neat et al., 2014). Although the Northwest component shows genetic homogeneity with the South, their connectivity is limited (Neat et al., 2014), and these units are often considered as separate populations. The focus of this study was primarily on exchange between the South and Viking units, therefore, the Northwest, Skagerrak, and Kattegat components were included solely as sink of larvae, not as sources.

Figure 1.

Spawning areas (numbered black polygons, Sundby et al., 2017) and population units (colour coded areas, Neat et al., 2014; ICES, 2015) included in the study. We only model cod eggs and larvae spawned within the South and Viking populations, but also estimate contribution to the other populations. Ocean regions outside the areas included in the study are hatched in red.

Due to both increasing temperatures and high fishing mortality, the North Sea cod stock declined substantially between the 1980s and early 2000s (Engelhard et al., 2014). Since Atlantic cod in the North Sea is close to the south-eastern border of its distributional range (Drinkwater, 2005), effects of increased ocean temperatures on larval survival, both directly and indirectly through, e.g. food availability, are likely to be strong (Beaugrand et al., 2003; Pörtner and Farrell, 2008; Akimova et al., 2016). To meet the management challenges of a changing North Sea cod population landscape, more research is therefore needed on both population connectivity and environmental effects on larval survival. Biophysical models will continue to be used for this purpose, and understanding the robustness of model results is therefore important.

Material and methods

To simulate the drift of cod eggs and larvae we coupled an ocean circulation model to a Lagrangian particle tracking model with an integrated cod egg and larvae individual-based model (IBM). We aimed to investigate the influence of (i) implementing realistic larval behaviour, (ii) improving ocean model resolution, and (iii) running simulations for different years on drift patterns and retention of larvae within population units.

Ocean models

We used the Regional Ocean Modelling System (Shchepetkin and McWilliams, 2005), a free-surface, terrain-following, primitive equations ocean model commonly used in biophysical modelling studies. We compared two ROMS reanalyses specifically designed for ocean regions adjacent to Norway, including the North Sea: SVIM, a relatively coarse resolution setup with 4 km horizontal resolution, 32 vertical layers, and output stored daily (Lien et al., 2013); and KINO, a finer resolution setup with 1.6 km horizontal resolution, 40 vertical layers, and output stored every 3 h (Sundby et al., 2017). To resolve light variation throughout the day, output from both reanalyses was interpolated to 1 h time steps in the particle tracking model (below). The SVIM reanalysis covers a wide region in the Northeast Atlantic and Arctic Ocean, while the KINO reanalysis covers a smaller area including the North Sea (Supplementary Appendix S1).

The SVIM reanalysis used the Simple Ocean Data Assimilation reanalysis version 2.1.6 (Carton et al., 2000; Carton and Giese, 2008) for initial and boundary values. Atmospheric forcing was taken from NORA10 (NOrwegian ReAnalysis, 6 h temporal resolution, 10 km horizontal resolution, Reistad et al., 2011), a dynamic downscaling based on the ECMWF (European Centre for Medium Range Weather Forecast) analysis. Shortwave and longwave radiation were analytically calculated internally. The internal time step was 150–180 s. The Generic Length Scale (GLS) scheme was used for parameterization of vertical turbulent mixing (Umlauf and Burchard, 2003). See Lien et al. (2013) for additional detail.

The KINO reanalysis used lateral boundary forcing from the GLORYS2V3 reanalysis (monthly average values, marine.copernicus.eu) consisting of a global, 1/4° Mercator grid, 75 vertical levels, 1 m top level, and 200 m bottom level (ORCA025). Atmospheric forcing was taken from the ERA-interim reanalysis (6 h temporal resolution, 0.25° horizontal resolution, Dee et al., 2011). The internal time step was 40 s. The Mellor Yamada closure scheme was used for vertical turbulent mixing (Mellor and Yamada, 1982). See Sundby et al. (2017) for additional detail.

Thus, in addition to spatial and temporal resolution, the reanalyses differ in the boundary conditions, atmospheric forcing, area coverage and vertical turbulent mixing scheme. However, since we focus on the North Sea away from the boundaries, we assume that ocean dynamics are mainly dominated by the internal ocean dynamics calculated by ROMS, and will therefore not differ systematically apart from due to resolution.

Simulating egg- and larval drift

To model the drift of planktonic eggs and larvae, we used the open source particle tracking framework OpenDrift (Dagestad et al., 2018, github.com/opendrift), which was coupled offline to the SVIM or KINO reanalyses. OpenDrift contains modules for estimating drift of planktonic eggs, and we additionally developed a cod egg and larvae IBM (github.com/trondkr/KINO-ROMS/tree/master/ICESJMS-2018), which was integrated as a module to OpenDrift. Eggs and larvae were advected horizontally using an Euler scheme which, given the small time step (1 h), showed negligible difference compared to more computationally expensive Runge-Kutta schemes (Supplementary Appendix S2). Horizontal diffusion was not included because it would introduce randomly driven differences between simulations not attributable to vertical behaviour, ocean model resolution or interannual variation.

The planktonic egg phase was parameterized for cod as a function of development time (D, days) dependent on ambient sea water temperature (T, °C) obtained from the ocean model reanalyses (Langangen et al., 2014 based on Ellertsen et al., 1987):

$$\text{ln}D=3.65–0.145\times T$$

(1)

After completing the egg stage, simulated individuals hatch into cod larvae. The cod larvae IBM contains modules for growth, vertical behaviour and mortality and was developed based on earlier modelling studies of larval cod (Kristiansen et al., 2009a, 2009c, 2014a). The growth rate of larvae (GR, percentage of larval weight day⁻¹) depends on larval weight (W, mg) and ambient temperature (T) (Folkvord, 2005), with initial weight set at 0.08 mg:

$$GR=1\mathrm{.0}8+1.79\times T–0.0\mathrm{74}\times T\times \text{ln\hspace{0.17em}}W–0.0\mathrm{965}\times T\times \text{ln}{W}^2+0.0112\times T\times \text{ln\hspace{0.17em}}{W}^3$$

(2)

Eggs and larvae are subject to vertical turbulent mixing parameterized from wind speed (Sundby, 1983), using a binned random walk scheme (Thygesen and Ådlandsvik, 2007; Dagestad et al., 2018). In addition, larvae exhibit vertical swimming behaviour in response to environmental conditions, representing the trade-off between feeding opportunity and predator avoidance. Specifically, larvae swim up if ambient light is decreasing and predation risk can be assumed to decrease, and down if light is increasing and predation risk can be assumed to increase. Swimming speed (SS, mm s⁻¹) is dependent on larval length (Peck et al., 2006):

$$SS=0.261\times (L{1.552}^{\times }{L}^{\mathrm{-0.08}})-(5.289/L)$$

(3)

and the fraction of each time step spent swimming is set to 15% (to represent pause-travel behaviour, Munk, 1995). Larval length (L, mm) is a function of weight (Folkvord, 2005):

$$L=e^{2.296+0.277\times \text{\hspace{0.17em}ln}W-0.005128\times \text{\hspace{0.17em}ln}{W}^2}$$

(4)

Eggs and larvae are subject to size-dependent mortality (m, daily). Similar to Akimova et al. (2016), we set egg mortality fixed at 0.2, and larval mortality decreasing with weight according to:

$$m=0.0\mathrm{6}\times W^{-0.4}$$

(5)

Instead of removing individuals, the survival-probability of each individual is updated throughout the simulation according to Equation (5). See Figure 2 for an overview of the IBM functions.

Figure 2.

Egg and larvae IBM functions. (a) Egg development time (D) as a function of temperature; (b) Daily growth rate (GR, contours) as a function of larval weight and temperature; (c) Swimming speed (SS) as a function of larval length; (d) Larval length (L) as a function of weight; (e) Mortality rate (m) for larvae as a function of larval weight; and (f) Depth distribution of a single, representative individual from the egg stage (dashed line) through the larval stage (solid line) to settlement.

We released particles representing cod eggs in the following main spawning areas: Dogger Bank, Dogger Bank Central, German Bight, Norwegian Trench, and Viking Bank, using recently updated spawning grounds for North Sea cod (Sundby et al., 2017) (Figure 1). The spawning period for North Sea cod lasts from around January to April (Sundby et al., 2017), likely happening earlier in the south than the north (Brander, 2005; ICES, 2017). However, data referred to in Heath et al. (2008) show little difference in peak spawning timing between central/northern areas (9–10th of March) and southern areas (3rd of March). We therefore released individual cod eggs from 15th of January to 15th of April at all spawning grounds, with the number of eggs released per day following a Gaussian curve peaking on 1st of March and the total number of eggs released per spawning ground per year summing to around 10 000.

At release, eggs are distributed uniformly within the spatial extent of the spawning grounds, at depths varying between 0 and 30 m at 10 m depth intervals. In runs without vertical movement, the vertical mixing and swimming behaviour schemes are deactivated, and eggs and larvae drift at the release depth. If an individual encounters a bottom depth shallower than its release depth, it continues to drift at the bottom depth.

We ran simulations using ocean model reanalyses for 2012 and 2013. For both years, the simulated drift of eggs and larvae was run until the 15th of August. We assumed that larvae could settle anywhere, and larvae that reached >49 mm during the simulation were considered as successfully settled (Bastrikin et al., 2014).

Analyses

We ran eight simulations (Table 1) to investigate the influence of (i) fixed depth vs. vertical movement of eggs and larvae, (ii) ocean model resolution (SVIM vs. KINO), and (iii) interannual variation in ocean dynamics (2012 vs. 2013). To compare the results, we created a 10 × 10 km binned model domain between −4 and 14°E and 53 and 62°N (approximately outlining the North Sea), and summed the number of settling larvae (>49 mm) within each bin. To test if results differed depending on the spatial and temporal scale of the analyses, we also compared results at hatching and metamorphosis (12 mm); and at coarser bin resolutions (50 and 300 km). The particle distributions for the five spawning grounds (Figure 1) were combined into a total cumulative particle density distribution (PDD). PDDs were weighted by the survival probability of the larvae (i.e. higher survival probability receiving higher weight).

Table 1.

Different model runs tested.

#	Model run	Year	Horizontal resolution (km)	Temporal resolution	Vertical movement
1	SVIM 2012 NV	2012	4	Daily	None
2	SVIM 2012 V	2012	4	Daily	Active
3	SVIM 2013 NV	2013	4	Daily	None
4	SVIM 2013 V	2013	4	Daily	Active
5	KINO 2012 NV	2012	1.6	3 h	None
6	KINO 2012 V	2012	1.6	3 h	Active
7	KINO 2013 NV	2013	1.6	3 h	None
8	KINO 2013 V	2013	1.6	3 h	Active

SVIM: Coarser ocean model reanalysis; KINO: Finer ocean model reanalysis; NV: No vertical movement; V: Vertical movement included. The name and number of model runs are used throughout the text.

We cross-compared PDDs using the fraction of unexplained variance (FUV) 1 – r², where r is the Pearson correlation coefficient of two vectorized PDDs (similar to Simons et al., 2013); and Fuzzy Kappa Index (FKI), a categorical map comparison technique recently applied for coupled biological-oceanographic systems (Rose et al., 2009; Stow et al., 2009). FKI is based on Kappa, a cell-by-cell comparison where numerical values are converted into categories (e.g. high, medium, and low number of particles). Kappa is computed based on a misclassification matrix, which sums the number of cells with category disagreement, and ranges from 1 (perfect agreement) to −1 (complete disagreement), where 0 indicates expected agreement between two uncorrelated maps. FKI extends this approach with fuzzy set theory, taking into consideration the neighbourhood of a cell (Hagen-Zanker et al., 2005; Hagen-Zanker, 2009). That is a cell is primarily defined by its category, but also by the categories in its neighbourhood. FKI also considers “fuzziness”, i.e. the distinction between categories can be gradual. Thus, compared to exact cell-by-cell comparisons as FUV, FKI also credits near cell-to-cell agreement, allowing distinguishing between small and large disagreement in position and category (van Vliet et al., 2013).

In the parameterization of FKI, the user defines the neighbourhood, categorical fuzzy sets and categories. Since no prior knowledge was available, we used unreasoned, sensible base values and ran sensitivity tests of parameters (Supplementary Appendix S3). We used the freely available software Map Comparison Kit to calculate the FKI (www.riks.nl/mck, Visser and De Nijs, 2006).

In addition to comparing PDDs from different model runs, we calculated connectivity matrices between population units, using the estimated geographical population extents from Neat et al. (2014) and ICES (2015) (Figure 1). Here, the Viking unit encompasses the spawning grounds Viking Bank and Norwegian Trench, while the spawning grounds Dogger Bank, Dogger Bank Central, and German Bight are included in the South unit. In addition to quantifying retention and transport within these two units, we differentiated export to Skagerrak, Kattegat, Northwest (east of Scotland) and outside the study domain (to the north and west).

Results

Patterns of modelled larvae distributions

Across all model runs, high numbers of larvae settle in southern and central parts of the North Sea, along the Scandinavian coastlines and, more dispersed, in northern parts of the basin (Figure 3). Also, considerable drift of larvae occurs into Skagerrak, but only a small amount of larvae arrives to Kattegat. Settlers in the northern North Sea are in 2012 concentrated around Shetland Islands and northern Scotland, and in 2013 more to the east in the Viking Bank area. In the southern North Sea, the bulk of settling larvae are concentrated around the central (Dogger Bank) area in 2012. In the SVIM runs, a high number of larvae also settle in the German Bight. In 2013, the SVIM runs result in lower numbers of settling larvae in the south compared to the KINO runs, where a high number of settlers are spread along the southern coast of the basin, with a peak in the German Bight. Survival probability is generally low, as eggs and larvae drift for a long period (130–213 days) before reaching settlement length. Still, since egg and larval duration decreases as temperatures increase in spring and summer, 98–100% of the larvae reach settlement before the 15th of August.

Figure 3.

Distribution of larvae at settlement in different model runs, pooled for all spawning grounds. Colour scale: Number of settling larvae within each 10 × 10 km cell, weighted by survival probability (i.e. particle density distributions, PDDs). Note that for illustrative purposes, the colour scale is square root transformed. NV: No vertical movement, V: Vertical movement included.

Quantitative comparison of scenarios

To quantitatively compare PDDs at settlement (Figure 3), we calculated the FUV and FKI between all model runs. Both metrics indicate that the largest difference between runs (i.e. highest FUV and lowest FKI) is driven by interannual variation in ocean dynamics (Tables 2 and 3). The most similar PDDs resulted from comparing runs with or without vertical movement, but for the same ocean model resolution and year. To further quantify the importance of vertical movement, ocean model resolution and interannual variation, we calculated the mean FUV and FKI value for all relevant comparisons, e.g. the effect of vertical movement was calculated as the mean of the comparisons for runs 1 and 2, 3 and 4, 5 and 6, and 7 and 8. The mean FUV and FKI values similarly indicate that the effect of interannual variation is highest and the effect of vertical movement lowest, while the effect of ocean model resolution is close to the effect of vertical movement (Tables 2 and 3, upper right). For the FKI, sensitivity analyses and analyses across levels of stage and scale (below and Supplementary Appendices S4 and S5) additionally show that the effect of ocean model resolution is always higher than vertical movement and often intermediate between interannual variation (high effect) and vertical movement (low effect). As expected, absolute FKI values are higher when larger neighbourhood is considered (Koch et al., 2015). Still, the relative patterns are maintained, and results are relatively insensitive to changes in fuzziness of category.

Table 2.

Fraction of unexplained variance (FUV) from comparisons of PDDs of settled larvae in different model runs.

Higher values and dark shading indicate dissimilarity, low values and light shading indicate similarity. Also shown (upper right): mean effect of vertical movement (VM), interannual variability (IV) and ocean model resolution (Re); and (lower right) sensitivity to vertical movement and interannual variability by ocean model resolution. NV: No vertical movement, V: With vertical movement.

Table 3.

Fuzzy Kappa Index (FKI) for comparisons of PDDs of settling larvae in different model runs.

See Table 2 for details. Note that values are inverted compared to FUV but shading scale is maintained: higher values and light shading indicate similarity, lower values and dark shading indicate dissimilarity.

Looking into detail, the FUV and FKI values indicate more different results when comparing runs with and without vertical movement for SVIM relative to KINO (Tables 2 and 3, lower right), suggesting that the coarser model is more sensitive to inclusion of vertical movement. This result is robust across FKI sensitivity tests and for different spatial and temporal scales of analyses (Supplementary Appendices S4 and S5). Conversely, sensitivity tests and comparison across scales show no clear difference between the two ocean models in sensitivity to interannual variation (Supplementary Appendices S4 and S5).

Patterns across temporal and spatial scales

The general patterns in FUV and FKI values are maintained at different temporal scales (calculating PDDs at hatching, metamorphosis, and settlement, respectively) and spatial scales of binning (10, 50, and 300 km, respectively) (Figure 4); with largest differences between years, followed by ocean model resolution, and lowest differences attributed to vertical movement (Table 4, See also Supplementary Appendix S5).

Table 4.

FKI comparisons of PDDs across spatial and temporal scales.

Values indicate the effects of vertical movement, ocean model resolution and interannual variability (rows, calculated by averaging FKI across 4 comparisons). Higher values and light shading indicate similarity (low effect), lower values and dark shading indicate dissimilarity (high effect). See Supplementary Appendix S5 for full results.

Figure 4.

Distribution of settled larvae across different temporal and spatial scales. The upper panels show the PDDs from run 1 (SVIM NV 2012, 10 × 10 km bins) at different temporal scales: (a) hatching, (b) metamorphosis, and (c) settlement. Spawning grounds are outlined in black. The lower panels show the PDDs for run 1 at settlement using different bin sizes: (d) 10 km, (e) 50 km, and (f) 300 km. Note that the colour scale differs between plots, and is square root transformed for illustrative purposes. “c” and “d” correspond to Figure 3a.

As expected, all variables become more important with longer model runs (i.e. at settlement relative to at hatching or metamorphosis). For example, while the influence of vertical movement at hatching and metamorphosis is limited, FKI values are substantially reduced at settlement, indicating a clear influence of vertical movement at this stage (albeit being less important than the two other factors). While effects of interannual variability or model resolution are not clearly related to spatial scales, the influence of vertical movement is reduced at larger bin size. Still, the FKI value is relatively low using 300 km bin cells at settlement, indicating that vertical movement affects settlement patterns at scales comparable to those of population dynamics.

Retention within population units

Across all model runs, retention at settlement is higher in the South than the Viking population unit (Table 5). The largest differences in retention rates in the South unit occur between years, with higher retention in 2013 than 2012 (see also Supplementary Appendix S6). A substantial fraction of larvae from South is transported to Skagerrak, and this export is higher in 2012 than 2013. For the Viking unit, a large fraction of larvae is in 2012 transported to the north-western North Sea, but in 2013, almost all export happens to the north outside the study area. Comparing results for different ocean model resolutions or with or without vertical movement show few clear trends, with the exception of retention within the Viking unit which tends to be higher at coarser ocean model resolution.

Table 5.

Population connectivity across model runs.

Numbers indicate the fraction of settling larvae originating from different population units (columns) that settle within different units (rows) (Figure 1). Spawning grounds Viking Bank and Norwegian Trench were combined into the Viking population unit and Dogger Bank, Dogger Bank Central and German Bight into the South unit. We also differentiated drift to Skagerrak, Kattegat, Northwest (East of Scotland), or other areas to the north or west of our study area (“Outside”). Dark shading indicates high retention within the population or high export to other units (high connectivity) while light shading indicates low retention/low connectivity.

Discussion

An increasing number of studies use individual-based biophysical models to explore dynamics of early life stages of fish and invertebrates, but sensitivity analyses are rarely performed (Peck and Hufnagl, 2012). One constrain might be computational cost. Our study aimed to identify where to most efficiently invest computational cost when modelling larval drift and population connectivity, using cod in the North Sea as a case study. Instead of comparing model results to observations, we cross-compared the output of different model runs to assess the sensitivity of the results to (i) inclusion of vertical movement, (ii) ocean model resolution, and (iii) interannual variation in ocean dynamics. While this does not imply which results are most realistic, we can assess the relative impact of each factor, and consider their importance depending on the temporal and spatial scale of the study.

We found that interannual variation had the largest influence on the results, confirming previous studies on the importance of interannual variation in ocean circulation for modelled fish larvae drift in the North Sea (Bartsch et al., 1989; Bolle et al., 2009; Dickey-Collas et al., 2009; Savina et al., 2010; Lacroix et al., 2013). Interannual variation in ocean transport has been attributed to variation in wind patterns (Bartsch et al., 1989; Bolle et al., 2009; Savina et al., 2010) and, ultimately, the North Atlantic Oscillation (NAO) index, with e.g. higher drift into Skagerrak (Jonsson et al., 2016) and to the western Norwegian coast and continental shelf (Huserbråten, 2017) during high-NAO conditions. We also found that retention within the South unit was lower in 2012, a year with positive anomalies of the NAO winter index (National Center for Atmospheric Research Staff, 2017), compared to 2013, a NAO-negative year. Previous modelling studies found that retention of cod larvae was higher for the South compared to the Viking unit (Heath et al., 2008; Huserbråten, 2017). Our results support higher retention rates in South, but imply that interannual variation in ocean dynamics plays a critical role in determining the degree of retention in both areas.

Second, ocean model resolution had a larger influence on the results than vertical movement. Based on these results, one should prioritise to run simulations for multiple years and use higher resolution ocean models if computational trade-offs are needed (as in e.g. Jonsson et al., 2016; Barth et al., 2017). However, depending on the model system, increasing resolution might be computationally costlier than implementing vertical behaviour. In our case, simulating larval drift from one spawning ground (∼10 000 particles, 211 days in 1 h time steps) on a laptop with 2.9 GHz dual-core Intel Core i5 processor took ∼4 h 15 min, which increased to ∼5 h 20 min by adding vertical movement or ∼10 h 30 min by using high-resolution KINO forcing instead of SVIM (without vertical movement). Running high-resolution ocean models to provide circulation data is also costly, e.g. running 1 month on a 512 CPU server took 2 h for SVIM but 3.5 days for KINO (Vidar Lien, pers. comm.).

Moreover, the importance of considering vertical movement depends on the spatial and temporal scale of the study. For example, despite being less important than the two other factors, the effect of vertical movement was non-negligible for small-scale patterns (distributions at 10 km scales). Moreover, at settlement after several months of drift, vertical movement also influenced large-scale patterns (distributions at 300 km scales). This is consistent with the large number of studies emphasizing the importance of including realistic vertical behaviour of pelagic larvae in biophysical models (e.g. Paris and Cowen, 2004; Cowen et al., 2006; Fox et al., 2006; Fiksen et al., 2007; Vikebø et al., 2007; Ospina-Alvarez et al., 2012; Drake et al., 2013).

The limited contribution of vertical movement compared to the two other factors is likely attributable to characteristics of the North Sea, a shallow shelf sea where strong tidal currents result in a well-mixed water column, notably in the southern region (Sundby et al., 2017). Other modelling studies have similarly found limited effects of vertical movement on drift of spring-spawned sole larvae (Savina et al., 2010) and winter-spawned herring larvae (Dickey-Collas et al., 2009) in the North Sea. Dickey-Collas et al. (2009) hypothesized that this was caused by the well-mixed and variable dynamics of the southern North Sea in winter and spring, and that the situation might differ for species spawning in fall when waters are more stratified. Importantly, while vertical behaviour may have a relatively low influence on larval drift patterns in shallow, well-mixed areas such as the North Sea, it will still be critical for food availability and larval survival (Fiksen et al., 2007; Kristiansen et al., 2009b).

While we are unaware of studies describing vertical movement of North Sea cod larvae, simulations using the same swimming speed function (Equation 3) corresponded well with data from other regions (Kristiansen et al., 2009a, 2014b). However, vertical behaviour was in these studies also driven by prey fields, generally confining larvae to the upper 40 m. Thus, our simulations with fixed depths < 30 m might more closely resemble larval depth ranges (despite ignoring individual movement), while our vertical movement scheme might overestimate this range.

Relative to previous efforts to compare drift model results using different ocean models, both SVIM and KINO can be considered as “high resolution”. Hufnagl et al. (2016) found significant inter-model variability, in some cases exceeding interannual variability, when comparing fish larvae dispersal using different ocean models, but models with horizontal resolution <10 km gave relatively similar results. Comparing a model with 0.08° horizontal grid and daily time step with coarser products (0.24° and 0.56° grid, 5 and 30 days), Putman and He (2013) found that coarser resolution generally resulted in higher offshore transport, likely by averaging out frontal zones between water masses that promote retention. Overall, we did not observe lower retention rates at coarser resolution, indicating that the reanalyses used capture relatively similarly features promoting retention.

Still, we observed differences between SVIM and KINO. With a horizontal resolution of 1.6 km vs. 4 km, KINO resolves smaller and/or more eddies than SVIM, which results in more horizontal spreading of the larvae. Temporal resolution may also cause differences, specifically, tidal forcing was included in both reanalyses (Lien et al., 2013; Sundby et al., 2017), but tides will only be resolved in the 3 h KINO output, not the daily averaged current fields from SVIM. Tides are likely important drivers of fish larvae transport in strongly tidally influenced regions such as the southern North Sea (Fox et al., 2006; Sundby et al., 2017). For example, it has been suggested that flatfish larvae perform tidally cued vertical swimming to increase retention in nursery areas (De Graaf et al., 2004), and resolving tidal forcing increased retention of modelled cod larvae in the Gulf of Maine (Huret et al., 2007). Consequently, it is generally advisable to use ocean models with sufficient temporal resolution to capture tidal variations when modelling larval drift in tidally influenced areas.

Differences in atmospheric forcing between the reanalyses could also influence drift patterns. In particular, since ambient ocean temperature drives development and growth rates of eggs and larvae in our simulations, which in turn determines larval swimming speed and mortality, consistent differences in temperature fields could substantially influence the results. However, we did not observe a consistent bias when comparing surface and depth-integrated temperatures between the reanalyses (Supplementary Appendix S7).

Finally, we found that with SVIM, results were more sensitive to inclusion of vertical movement. Monthly averaged current fields for July (Supplementary Appendix S8) show that general patterns are comparable in the two reanalyses, and the water column is as expected relatively well mixed. Exceptions occur along the coastal currents north of Denmark and along the southern Norwegian Coast, where surface currents are much stronger than the vertically integrated current fields. In particular along the Norwegian coast, this vertical gradient is stronger in SVIM than KINO. As larvae drifting passively in the upper layer are more likely to be transported with surface currents, this may explain some of the observed differences. Differences in vertical gradients might be driven by the turbulent mixing schemes, i.e. the Mellor Yamada scheme used in KINO potentially smooths the vertical layers more than the GLS scheme used in SVIM.

The connectivity rates estimated here appear sufficient to homogenize some of the North Sea cod populations. High connectivity from the North Sea to Skagerrak is supported by genetic data (André et al., 2016), but despite moderate estimated transport levels from South to Viking, these populations are genetically distinct (Heath et al., 2014). Importantly, population isolation is not solely driven by larval retention as assessed in this study, but also by natal homing (Neat et al., 2014; André et al., 2016), and, potentially, by selection for locally adapted phenotypes (Barth et al., 2017). Moreover, realized larval retention and connectivity will not only depend on drift, but on available settlement habitat and on post-settlement survival, driven by factors such as predation, food availability and competition with other settlers (Heath et al., 2014). Taking these factors into account demand more complex biophysical models and fine-scale observation data, and was beyond the scope of this study. However, it should be considered in future studies aiming to realistically quantify larval connectivity and retention.

Observed changes in the distribution of North Sea cod have been attributed to local overexploitation and increasing temperature, leading to a decrease of the South unit (Engelhard et al., 2014). But while the South and Viking units show diverging population dynamics and display genetic and ecological differences, the role of metapopulation connectivity is only partially understood. For example, little is known about long-term variability in larval retention within population units. Similarly, a large body of literature exists on metapopulation connectivity, in particular from coral reefs, but relatively few studies have investigated multi-year dynamics (reviewed in Cowen and Sponaugle, 2009; Jones et al., 2009). Exceptions include studies finding largely consistent retention and connectivity patterns using parentage analysis (Saenz-Agudelo et al., 2012) or geochemical tags (Carson et al., 2010), or significant interannual variation using genetic assignment tests (Hogan et al., 2012). Based on the present and previous biophysical modelling studies, it is critical that future studies take interannual variation into consideration, in particular when considering the potential implications of these studies for conservation and management (Botsford et al., 2009; Heath et al., 2014).

Acknowledgements

K.Ø.K. was supported by the Research Council of Norway (RCN) project SUSTAIN (244647/E10) and the WHOI John Steele Postdoctoral Award. G.R. was supported by the Norden Top-level Research Initiative sub-programme “Effect Studies and Adaptation to Climate Change” through the Nordic Centre for Research on Marine Ecosystems and Resources under Climate Change (NorMER). Ø.L. was supported by the RCN project OILCOM (255487). Simulations were performed on resources provided by UNINETT Sigma2—the National Infrastructure for High Performance Computing and Data Storage in Norway. We thank two anonymous reviewers for their useful comments on the manuscript.

Author notes

The Kristina Øie Kvile and Giovanni Romagnoni authors contributed equally to the article.