Two Decades of Match-Mismatch in Northeast Arctic Cod – Feeding Conditions and Survival

The successful recruitment of Northeast Arctic (NEA) cod is thought to depend on sufficient and suitable prey for the newly hatched larvae, in particular the nauplii stages of the lipid-rich calanoid copepod species Calanus finmarchicus. The role of spatial and temporal variations in prey availability in combination with temperature and other factors in influencing growth and survival of cod larvae is, however, incompletely understood. By combining an individual based model for NEA cod larvae at the Norwegian coast with a high-resolution ocean model and a nutrient-phytoplankton-zooplankton-detritus model providing 18 years of daily environmental conditions and prey availability we assessed larval growth and survival until they settle in their feeding habitat in the Barents Sea in early fall. We find on average a two-week delay from the peak timing of first-feeding cod larvae to the peak in prey availability. In warm years, more larvae experience food limitation than in normal years. The positive effects of high temperature on growth, survival and ultimately recruitment are nonetheless larger than the negative effects of food limitation. Food limitation mainly affects larvae spawned in southern areas or late in the spawning season as these larvae experience the highest temperatures and have the highest energy requirements. Our findings highlight the spatial and temporal differences in mechanisms that regulate growth and survival of early life stages of NEA cod and suggest that spatially resolved data may be essential for understanding match-mismatch dynamics.


INTRODUCTION
Variability in the number of individuals in adult fish populations results mainly from the loss of individuals through natural and fishing mortality, and recruitment of individuals through the reproduction process (Payne et al., 2009). The survival until recruitment to the harvestable stock is a complex process that encompasses all conditions experienced by the fish including the variability in physical and biological conditions (Hjort, 1914). Numerous hypotheses have been developed to explain the observed variability in year-class success, mostly related to the early life stages of fish development. The "Critical Period" hypothesis explains year-class survival variability through the success or failure in finding enough suitable prey after the yolk-sac has been exhausted. Meanwhile, the "Aberrant Drift" hypothesis states that larvae survival is dependent on the dispersal to favorable areas, i.e., with enough prey and suitable abiotic conditions (Hjort, 1914;Houde, 2008). The "Match-mismatch" hypothesis explains the year-class survival variability through phenology and extends the critical period from the first-feeding stage until after metamorphosis (Cushing, 1969(Cushing, , 1990. This hypothesis highlights mechanisms that connect the temporal overlap between prey and predators as a key for the survival of the predators, assuming spatial overlap is present. Mismatch with prey results in slower growth and development and may lead to enhanced mortality because the larval fish are more susceptible to predation (Hjort, 1914;Durant et al., 2007;Ferreira et al., 2020). The three above mentioned hypotheses account for many underlying processes of various spatial and temporal resolutions that affect survival of the early life stages of fish, including starvation, transport to and from favorable habitats, and metabolic rates due to varying ambient conditions. Moreover, the relative importance of these processes to recruitment varies in space and time (Gallego et al., 2012;Hare, 2014).
Natural variability and ongoing climate change significantly influence the spatial and temporal distributions of predators and prey and may alter phenology and trophic interactions (Durant et al., 2007). It can also modify the costs and benefits associated with different spawning times and areas, and cause constraints at spawning grounds and shifts in the feeding areas, potentially through climate-driven changes in prey distribution (Fossheim et al., 2015). Changes in temperature affect the habitat selection of the spawners (Sundby and Nakken, 2008;Höffle et al., 2014;Sandø et al., 2020). Offspring transport may also be affected by changing ocean currents and cross-shelf exchange due to changes in ambient temperature, wind patterns and stratification that affects vertical egg distribution and spring bloom dynamics through nutrient availability (Vikebø et al., 2010(Vikebø et al., , 2021Skagseth et al., 2015;Strand et al., 2017).
Demographic structure and biological condition of the spawners may also influence the survival and distribution of the eggs and larvae of fish. The age and size composition of the spawning stock can influence spawning location, time, and duration (Sundby and Nakken, 2008;Opdal and Jørgensen, 2015;Langangen et al., 2019). For example, when there are older and bigger fish in the spawning stock of Northeast Arctic (NEA) cod, there is higher egg abundance and a spatially more widespread egg distribution (Stige et al., 2017). Furthermore, maternal characteristics affect egg quality, which in turn can impact larval susceptibility to starvation (Marteinsdottir and Steinarsson, 1998).
The diet of first-feeding larvae of NEA cod is mainly composed of copepod nauplii and for older larvae copepodites of Calanus finmarchicus (Dalpadado et al., 2009;Ottersen et al., 2014). Larvae are restricted to forage on prey smaller than their gape size, so it is important that the larvae are in a favorable place and time regarding prey concentrations of suitable sizes. The development of their main prey item C. finmarchicus is dependent on the onset of the seasonal spring bloom, which varies considerably (weeks) between years (Ellertsen et al., 1987;Vikebø et al., 2019;Dalpadado et al., 2020). Observations of prey availability are limited in space and time and can only provide a snapshot of the processes occurring in the ecosystem, but recently there are indications of a shift towards lower abundance, especially of C. finmarchicus, at the Norwegian continental shelf (Toresen et al., 2019).
Species-environment interactions are non-linear, and interannual dynamics cannot be reliably inferred from spatially averaged data (Ciannelli et al., 2007). On the other hand, individual based models (IBMs) coupled with ocean hydrodynamics models and nutrient, primary production, zooplankton, and detritus (NPZD) models represent complementary information about availability of suitable and sufficient prey items. The high spatial and temporal resolution of such models allow studying how local-scale processes translate into population-scale dynamics. In this study, we quantify availability of suitable prey items by combining the results from an IBM of NEA cod early life stages with an ocean hydrodynamics model and a NPZD model across the period 2000-2017. This approach enables quantification of intra-and interannual variations in match-mismatch and survival of NEA cod larvae through considerations of individual prey catches along individual drift trajectories for integration across the entire population (Kristiansen et al., 2009;Daewel et al., 2011;Strand et al., 2017). The analyses of the simulation results also enable investigation of the underlying spatio-temporal mechanisms regulating dispersal, feeding and development of individual NEA cod spawned along the Norwegian coast. We specifically address the following questions: (I) Does prey availability (timing, size, and abundance) affect individual larval growth, spatially, seasonally and interannually? (II) How does this effect on growth impact survival through early life stages until recruitment? (III) What are the key mechanisms and their relative importance that control growth and survival until settlement at favorable nursery grounds?

MATERIALS AND METHODS
Individual NEA cod larval growth was investigated using offline coupled models. Overall, the model results were calculated by using an ocean hydrodynamics model, a nutrientphytoplankton-zooplankton-detritus model including an individual based model for Calanus finmarchicus, and an individual based model for NEA cod eggs and larvae. The simulations were run from 2000 to 2017, for the period between the 15th of January to the 30th of September each year. The results from the simulations were used to analyse growth and survival of the individuals from the egg stage until the larvae reaches the Barents Sea. We also account for size-dependent survival through the first winter after settlement in the Barents Sea.

Model Species
The Northeast Arctic Cod The Northeast Arctic (NEA) stock of Atlantic cod (Gadus morhua) is one of the world's largest and most studied cod stocks due to its commercial importance and of high commercial and ecological importance (Yaragina et al., 2011;Ottersen et al., 2014). NEA cod inhabits the Barents Sea, a shallow shelf sea and a transition zone for warm and saline waters between the Atlantic Ocean and the Arctic. Mature NEA cod individuals migrate long distances during winter-spring from their feeding habitats in the Barents Sea towards the south to spawn along the Norwegian coast (Figure 1), from Finnmark (71 • N) to Møre (63 • N), with the highest spawning activity in the Lofoten archipelago (69 • N) Sundby and Nakken, 2008). The main spawning season lasts from early March until the end of April with a peak in early April (Ottersen et al., 2014). After spawning, the offspring drift northeastward in the upper layers of the water column and are transported by the Norwegian Coastal current (NCC) and the Norwegian Atlantic current (NAC) from the spawning grounds back to the Barents Sea (Vikebø et al., 2005). Egg development until hatching is highly dependent on ambient temperature, but it usually lasts from 2 to 5 weeks (Ellertsen et al., 1989;Yaragina et al., 2011). The NEA cod larvae yolk sac provides nutrition for about 5-7 days, whereafter the larvae switch to exogenous feeding about 3-6 weeks after spawning (on average late April) (Ottersen et al., 2014). By August -September successful offspring have reached the nursery grounds in the Barents Sea and are denoted 0-group. During the autumn (September -October) the 0-group descend towards the bottom and switch from the pelagic to the demersal phase in the Barents Sea. Surviving individuals are recruited to the harvested stock when they reach 3 years age.
important link between the phytoplankton and higher trophic levels (Melle et al., 2004), e.g., larvae of NEA cod that feed upon C. finmarchicus nauplii and copepodites (Ellertsen et al., 1989). The C. life cycle consists of an egg stage, 6 nauplii stages (denoted N1 to N6), 5 copepodites stages and the adult stage (male or female). C. finmarchicus usually has a one-year life cycle, but this species can develop more than one generation in warmer regions, as in the southernmost parts of the Norwegian Sea . The reproduction of C. finmarchicus is linked to the seasonal development of the phytoplankton (Melle and Skjoldal, 1998;Stenevik et al., 2007). C. finmarchicus release their eggs in the upper layers of the water column during spring and develop with the spring bloom to reach the overwintering stage (mainly copepodites stage five) in mid-late summer. Overwintering occurs in deep (>500 m) cold water and lasts until late winter, when the overwintered individuals ascend to surface waters to reproduce (Skjoldal et al., 1992).

Numerical Model Components
The Hydrodynamical Model Component Ocean physics were modeled using the Regional Ocean Modelling System (ROMS). ROMS is a three-dimensional free surface, baroclinic, primitive equations ocean model with topography following s-coordinates in the vertical (Shchepetkin and McWilliams, 2005). In the present study, the Nordic Seas numerical ocean model hindcast (SVIM) archive with 4 km horizontal resolution and 32 sigma coordinates in the vertical was used. It covers the Nordic, Barents and Kara seas and parts of the Arctic Ocean (for further information we refer to Lien et al., 2013Lien et al., , 2016.

The NORWegian ECOlogical Model (NORWECOM)
NORWECOM is a coupled physical, chemical, biological model system (Skogen et al., 1995;Skogen and Søiland, 1998) applied to study primary production, nutrient budgets and dispersion of particles such as fish larvae and pollution. In the present study it was run in offline mode using physics (light, wind, hydrography, horizontal and vertical movement of the water masses) from the SVIM archive. The NORWECOM model includes an ordinary NPZD model with processes such as primary production, respiration, algal death, remineralization of inorganic nutrients from dead organic matter, self-shading, turbidity, sedimentation, resuspension, sediment burial and denitrification. In addition, the model has an IBM for C. finmarchicus (Hjøllo et al., 2012).

The Calanus finmarchicus Individual Based Model
The three-dimensional IBM for C. finmarchicus includes feeding, growth, mortality, movement, reproduction, and adaptive traits that control the interaction with the environment. The entire life cycle of C. finmarchicus is included in the model, resolving the 13 life stages and their features, also including a genetic model algorithm which allows the model to evolve individual traits. Up until the third nauplius stage the development and stage longevity is based on temperature. From the third nauplius stage (N3) C. finmarchicus starts feeding and is 2-way coupled to the NPZD model. Growth is a function of phytoplankton density, temperature and size based on a bioenergetics model. A more detailed description of the model can be found in Hjøllo et al. (2012) and Huse et al. (2018).

The NEA Cod Egg and Larvae Individual Based Model
The IBM for NEA cod eggs and larvae is a Lagrangian dispersal model, and it includes dynamical vertical positioning of eggs dependent on buoyancy and turbulence, prey-and temperature dependent larval growth, predator-prey encounter rates and capture, larvae ingestion rate and metabolism, and larval vertical migration (Vikebø et al., 2005(Vikebø et al., , 2007Kristiansen et al., 2007).
The NEA cod IBM was run offline based on daily ROMS currents and hydrography, and NORWECOM 2D fields of the C. finmarchicus from the Lagrangian IBM in addition to 3D Eularian fields of mesozooplankton for the period January 15th until September 30th, 2000-2017.
Individual particles were initiated as eggs at 10 known spawning sites along the Norwegian coast (Figure 1 -colored circles; Sundby and Nakken, 2008). While the spawning season lasts from early March until the end of April, we initiated eggs from February 15th to May 15th to enable consideration of conditions also outside the main spawning period. For some of the analysis model results were adjusted according to the observed spawning intensities in space and time (Strand et al., 2017). Specifically, the relative importance of each individual particle was weighted by a gaussian function of spawning date with a peak on April 1st and a standard deviation of 14 days, and by year-and spawning site-specific weights that scaled the spatial relative importance of spawning sites according to egg surveys and observations of spawning adults (following Strand et al., 2017).
A total of one hundred eggs were released at each of the ten spawning sites (SS) per spawning date (SD), amounting to 90 000 particles per year. Particles were released initially as eggs at about 30 m depth. Individual particle position was updated every 30 min according to a 4 th order Runge-Kutta advection scheme and vertical egg positioning in the water column was based on the density difference between the eggs and the ambient water and turbulence (Thygesen and Ådlandsvik, 2007). Eggs were initiated with a buoyancy of 31.25 and a population standard deviation of 0.69 (Sundby, 1997). Individual hatching time followed Geffen et al. (2006)-see Supplementary Information (SI) for details.
After hatching, eggs continued as larvae with an initial weight of 45 µg and a length of approximately 4 mm (Folkvord, 2005). Larvae performed diel vertical migration at a third of a body length per second, ascending during night and descending during day (see Vikebø et al. (2007) for more details). Light decays exponentially with depth and depends on surface irradiance and the attenuation coefficient (see SI for details). Hence, each egg and larva experienced unique abiotic and biotic conditions through the period of pelagic drift.
After hatching, the larvae immediately started feeding in our model, and thus they did not go through the stage of yolk-sac larvae. It is estimated that yolk-sac cod larvae have a duration of 5-7 days in nature. Within this period, the larvae were feeding mostly on nauplii N2, and included nauplii N3 about 10-15 days after hatching depending on temperatures in the range 4-6 • C (Supplementary Figure 1).
Larval growth in the IBM was potentially food limited, depending on available prey concentration, the metabolic costs and egestion/excretion losses. Larvae fed on mesozooplankton from the NPZD model and from the C. finmarchicus IBM, which together provided prey field concentrations per stage [ind/m 2 ]. Unlike in the C. finmarchicus IBM, mesozooplankton in the NPZD model is not stage resolved. They were therefore converted to C. finmarchicus equivalents assuming a similar relative stage distribution (see SI for equations and details).
Cod larvae may feed on prey stages limited by their respective gape size, here 8 % of the larval length according to Otterå and Folkvord (1993). Feeding was parameterized according to Fiksen and MacKenzie (2002), already included in an IBM and tested for sensitivity in a numerical study by Kristiansen et al. (2007) -see SI for details. Stomach content [µg] was assumed to reach a maximum of 6 % of the body weight of the larva according to Lough et al. (2005).
Larval growth potential (i.e., growth assuming unlimited food) was limited by temperature according to Folkvord (2005) and may only be reached if sufficient prey is ingested (Supplementary Figure 1) -see SI for details.
The higher the ambient temperature, the higher the metabolic requirements. Consequently, a higher quantity of food is needed. If the food requirement is met the larvae may grow faster and prey composition more quickly diversifies (Supplementary Figure 1). Note that larval weight may decrease while larval length may not. This was also considered in the model.

Data Analyses
Annual averages of individual larval net distance from SS and ambient temperatures 30 days post-spawning were compiled in a scatter plot to indicate inter-annual varying conditions for NEA cod early life stages. In addition, we extracted annual average modeled temperatures during March-August in the upper 50 m of the water column in a polygon (Figure 1, black contour) containing 95 % of 15 days post-hatch (dph) larvae across all years. This calculation allowed inter-annual comparison of environmental conditions independent of egg and larva dispersal trajectories. The polygon extracted mean temperature was used to separate warm and normal years using the median value (Table 1). Here, we used the nomenclature normal years instead of cold years because the temperature anomalies (not shown) for the period from 1921 to 2018 based on the Kola section showed there are no actual cold years in our study period 2000-2017 (Yashayaev and Seidov, 2015).

Survival
Cumulative survival until the end of the simulation (September 30th) was calculated as a product of several factors. In short, survival at the end of the experiment was the product of survival through hatching, not starving (nonSt), reaching the Barents Sea (inBS), stage-dependent mortality, and size-dependent mortality.
To be considered alive the eggs must have hatched before the end of the simulation. Also, we required the larvae to reach the Barents Sea, the NEA cod nursery grounds, implemented as being east of 19 • E. Furthermore, we required that the individuals did not starve at any time throughout the drift period. The median temperature value, 6.27 • C, was used to separate normal and warm years.
Starvation was defined as weighing 20% less than its expected weight (i.e., only temperature dependent weight -potential weight) given its length (see Folkvord (2005) for function relating length and weight). As larvae cannot become shorter, this may occur if larvae are unable to ingest prey that at least match metabolic costs and assimilation losses.
We also calculated stage-specific survival for each individual that accounted for baseline natural mortality from, e.g., predation. Instantaneous natural mortality rates (M) for each stage were: for eggs (M egg ) 0.169 d −1 , for larvae (M lar ) 0.075 d −1 and for early juveniles (M eju ) 0.023 d −1 (Langangen et al., 2014b;Bogstad et al., 2016). Individuals were larvae if they were within 4.2 and 24 mm in length and early juveniles within 25 -99 mm (Herbing et al., 1996;Folkvord, 1997;Hamre, 2006). Cumulative baseline natural mortality as a function of time is therefore; where the different M's are daily instantaneous stage-specific mortality rates, N t is the abundance at time t and N 0 is the abundance at time 0 and t egg , t lar and t ejuv are the number of days spent in the different stages.
We also accounted for size-dependent survival through the first winter. Specifically, we calculated a proxy for size-corrected abundance at age 1 (N 1age ) based on the abundance (N 0age ) and mean length (L 0age ) at the end of the simulation. The sizecorrected abundance was calculated for each SS and SD using the following equation: N 1age = N 0age * exp(c 1 * [ln(L 0age )ln(L mean )]). Here, the coefficient c 1 = 9.78 was estimated by Stige et al. (2019, Table S6) based on year-class survival of NEA cod through the first winter. The scaling factor L mean was the mean length of all individuals through all the years in the simulation.

Food Sufficiency Index (FSI)
In addition, we quantified a food sufficiency index (FSI) to find which individuals might have suffered from food limitation but are not food limited enough to be considered starved. FSI reports the fraction of individuals that grew at rates above 90 % of the food-unlimited temperature dependent growth potential for more than 90 % of the time until the end of the simulation. The FSI was included as a more sensitive measure than starvation of food-limitations on larval growth through the planktonic stages.

Drift and Development of Egg and Yolk-Sac Larvae
Egg stage duration until hatching varies with ambient temperatures, which in turn depends on SS, SD, year, and individual drift trajectory. Individuals spawned late and to the south generally experience higher temperatures (∼ 7.0 ± 0.3 • C) than those spawned to the north (∼ 4.8 ± 0.2 • C) (Figure 2A), while inter-annual variability ( Figure 2B) is highest early in the season likely because of variation in the onset of seasonal heating. Individuals originating from SS 4 (inside the Lofoten archipelago) experience anomalously low ambient temperatures (∼ 4.9 ± 0.4 • C) ( Figure 2C) compared with other SS at similar latitudes, likely because of less influence of relatively warm Atlantic Water, retention, and runoff. Warm years generally result in less variable ambient temperatures between SSs, with 2012 as a good example.
Duration of the egg stage affects both the time and location in which the larvae start feeding. With average egg stage duration of 18.7 ± 5.3 (s.d.) days in our model and yolk sac absorption (not included in the model) of another 5-7 days the first feeding of larvae typically occurs within the first 30 days after spawning. Low ambient temperatures (< 4.5 • C) experienced by the offspring correspond to long distance traveled, while elevated ambient temperatures (> 4.5 • C) may correspond to a variety of mean distances from SS (Figure 3). The association between low ambient temperature and long distance traveled may be linked to a general trend of decreasing temperature with increasing distance from the spawning grounds. The absence of clear dispersal patterns for the individuals experiencing higher ambient temperatures is likely due to more variable circulation features. We believe that due to a constant NCC, low ambient temperatures result in a prolonged drift period until 30 days after spawning and a longer distance traveled from the SS. In warmer years (annual average polygon temperatures) the individuals generally experience higher ambient (5.2 ± 0.3 • C) temperatures along their drift routes than in normal years (4.6 ± 0.1 • C). In 2004, the individuals on average travelled the farthest

Feeding Conditions for Larvae
In our simulation results, the peak timing in the abundance of larvae at 15 dph in the polygon occurred on average about 15 days before the peak timing of prey of suitable size (Figure 4). The mean values were calculated for the main distribution area of cod larvae (the polygon in Figure 1) for warm (red) and normal (blue) years. The shading for larval and prey abundances shows the interannual variation (± 1 standard deviation) for warm or normal years. The lines have been smoothed by a moving mean of 5 days to emphasize trends. Abundance of 15 dph larvae appear on average 8 days later in normal than warm years (Figure 4, solid lines). Meanwhile, prey abundance of suitable sizes peaks almost at the same time in warm and normal years (Figure 4, dashed lines), but the abundance of prey is on average 31% lower in normal years than in warm years (non-significantly lower, t = 1.17, p = 0.28). C. finmarchicus nauplii at stage N3 constitutes a major part of the food of the correct size for cod larvae at this age ( Supplementary  Figure 3), and the seasonal curve for N3 abundance resembles the seasonal curve in prey of suitable size (Figure 4).
Match between larvae and suitable prey is higher for warm years than normal years (Figure 4 -area below intersection of solid and broken red lines vs. blue lines; for interannual variability see Supplementary Figure 3). Both normal and warm years (e.g., 2003 and 2002) may have higher abundances of N3 than the mean, and curiously, the year 2016 (warm) presented a delayed peak of N3 as compared to other years (Supplementary Figure 2B). Elevated abundances of C. finmarchicus progress northwards along the coast with time before becoming limited to only the southern part of the model domain from early June and onwards (Supplementary Figures 2A,C). In higher latitudes it is possible to observe two peaks of elevated C. finmarchicus N3 concentrations ( Supplementary Figure 2A -area I and II), whereas in the southernmost area III, C. finmarchicus N3 is available throughout the entire simulated months.

Mechanisms of Survival of the Early Juveniles
The FSI shows that on average early spawned larvae experience sufficient and suitable prey conditions to support growth according to their temperature potential ( Figure 5A). Later in the year, as temperature and thereby prey requirements increase, larvae experience growth limited by food availability, and this is particularly strong for late spawning far south. FSI is more variable for the individuals spawned far south and early in the spawning season, and at all spawning sites for those spawned around the main spawning period from late March to early April ( Figure 5B).
Mean ambient temperature along egg and larval drift routes and the FSI have a negative relationship ( Figure 6A). Consistently, most warm years, according to mean polygon temperature from March -August (Table 1), tend to have a lower FSI, while normal years tend to have a higher FSI ( Figure 6A). Note, however, that 2012 and 2017 are warm years, with median ambient temperatures along drift routes and high FSI. Larval length (and weight -not shown) against FSI ( Figure 6B) displays similar patterns as for ambient temperature against FSI (Figure 6A). Larval length (and weight -not shown) is positively correlated with ambient temperature (Figure 6C).  (Figure 7 -dark green line with diamonds). If also accounting for size-dependent survival through the first winter from age 0 to age 1 , survival decreases even further still, resulting in particularly low survival in the same years as before (Figure 7 green petrol line with filled diamonds). The calculated survival to age 1 and the field data for the 0-age group (age-0 fish per spawning stock biomass) have a non-significant weak negative correlation (r = −0.12, p = 0.64). However, although also non-significant, calculated survival to age one is weakly positively correlated with recruitment success at age three (recruits per spawning stock biomass, from assessment reports of the Arctic Fisheries Working Group (ICES, 2020), r = 0.29, p = 0.25). The recruitment at age 3 showed a long-term decline through the study period, which was not seen in the age-0 data or in the modeled survival (Figure 7). Short-term variability in survival to age-1 and recruitment, represented by year-to-year differences, were also non-significantly weakly positively correlated (r = 0.07, p = 0.78).
The latitudes of NEA cod larvae at 15 dph are significantly positively correlated with starvation (r = 0.51, p = 0.03) and with 0-age length (r = 0.52, p = 0.03). The strong correlation between larval latitude and longitude is because of the northeasterly orientation of the coast between spawning and nursery grounds.
The day of the year with the highest larvae abundance naturally correlates positively with their mean latitude (r = 0.83, p < 0.001) and longitude 15 dph (r = 0.71, p < 0.001) because FIGURE 5 | Fraction of individuals that grow 90% or more than their temperature-dependent potential if the prey is unlimited for more than 90% of the time since hatching until the 30th of September (food sufficiency index), per spawning site and spawning date.  a prolonged stage duration allows a longer drift period until the date of 15 dph. The time difference between peak abundance of 15 dph cod larvae and C. finmarchicus N3 is positively correlated with northward drift and day of peak abundance C. finmarchicus N3. Peak abundance of 15 dph larvae is on average 15 days prior to the peak of C. finmarchicus N3 and a delayed peak indicates a prolonged egg stage duration and therefore prolonged dispersal time to reach 15 dph. Latitude of the highest concentration of N3 has a positive correlation with cod survival to age 1 (r = 0.5, p = 0.035).

Key Findings and Implications
The results presented here shed light on the processes and mechanisms that affect larval growth, distribution and survival from NEA cod egg and larval stages until arrival of juveniles in their nursery grounds in the Barents Sea and through the first winter by considering previously established empiricalbased size-dependent survival prospects. Using a coupled ocean hydrodynamical model, a NPZD model including a C. finmarchicus IBM and an IBM for pelagic stages of NEA cod, it was possible to analyse the importance of abiotic conditions and food availability for growth and survival through all early life stages until settlement of juveniles at nursery grounds in the Barents Sea for 18 consecutive years. Observationalbased abundance estimates were also added at multiple stages and statistical correlations were performed with the modeled larval abundances. In comparison a similar study focused on hydrodynamical processes along the coast that affect spatialtemporal overlap between NEA cod larvae and prey for a few selected years (Vikebø et al., 2021). On the whole, the two studies provide complimentary information on key physical processes along the Norwegian coast that are important to NEA cod larvae feeding and to what extent combined biophysical models may predict the population size at different development stages.
We found that the time of peak abundance of 15 dph larvae is around two weeks before the peak timing of prey, with considerable interannual variability both in timing and amplitude. Because of this mismatch in peak timing, the NEA cod larvae seem to depend more strongly on the time of the increase in prey concentration than on the time of peak or decrease in prey. A delay between cod larvae and their prey is consistent with observations from Vestfjorden in the Lofoten archipelago, where an early peak of C. finmarchicus copepodite stage C1 in warm years tended to coincide with large mean size of NEA cod larvae around 15 dph, presumably because of a good temporal match between predators and prey (Ellertsen et al., 1987(Ellertsen et al., , 1989Cushing, 1990).
Most larvae grow at or close to their temperature potential except for late-spawned individuals. Late spawned cod eggs hatch into first-feeding larvae with better temporal match with suitable prey, at the same time as seasonal heating allows for a higher growth potential and thereby elevated feeding requirements to support growth and metabolic costs. In fact, our model results show that fast-growing larvae experience more often insufficient FIGURE 7 | Cumulative modeled survival until the 30th of September only including fixed stage-specific mortality rates (yellow line with upside-down triangles), also adding requirement that individuals need to reach the Barents Sea (green with triangles), also requiring that individuals do not starve (dark green with open diamond) and including size-dependent survival until age 1 (blue line with filled diamonds). For comparison we show observed 0-age group abundance divided by spawning stock biomass (purple line with filled squares) and abundance of 3-year olds from Virtual Population Analysis back-calculated to spawning year divided by spawning stock biomass (SSB -dark blue line with filled circles). Units for the latter two variables are thousands per tones, data from the ICES (2020).
prey conditions to match their growth potential depending on their size and ambient temperature. March -August mean temperature during the early life stages of NEA cod in the core area of appearance has positive correlations with the length and weight at 0-age (r = 0.80, p < 0.001; r = 0.88, p < 0.001, respectively) and a negative correlation with the FSI (r = −0.62, p = 0.005). This is also the reason why they are more prone to experience starvation as they have elevated metabolic costs and therefore require an elevated steady supply of suitable prey which in our model is not always available. However, the reward of finding sufficient food is rapid growth and high cumulative survival.
First-feeding and young larvae depend on a narrow selection of C. finmarchicus stages. The results here show that the abundances of prey for these stages are still significantly affecting NEA cod weight and survival at the 0-age group and 1year old individuals. Our results further suggest that this association at the population level is mainly driven by the food conditions experienced by a minority of individuals that have high growth potential because of high temperatures, as these are the ones that are potentially food-limited. Despite having sufficient food, the individuals that originate from early spawning and northern spawning grounds develop slowly and have low survival because of low ambient temperatures. Individuals that originate from late spawning and southern spawning grounds have higher growth potential because of higher temperature, but also higher risks of food limitation and of not reaching the Barents Sea nursery grounds. An implication of this spatial heterogeneity in which mechanisms regulate growth and survival is that the dynamics are not necessarily well represented by spatially averaged data (as pointed out by, e.g., Ciannelli et al., 2008). For example, average food and temperature conditions may not be the most relevant for yearclass survival, but rather food levels in the water masses with high temperature.
Dispersal to the north is necessary to reach the NEA cod nursery grounds but also impose an enhanced risk of mismatch because suitable C. finmarchicus nauplii abundance is only present at the shelf towards the entrance to the Barents Sea during a much more limited period than farther south. This is likely also the explanation why we find a positive correlation between larval latitude and starvation. During spring (April -May) C. finmarchicus N3 individuals are abundant both inside the polygon encompassing 95 % of 15 dph cod larvae and even farther north. However, toward summer (June-July) C. finmarchicus N3 are almost absent in the polygon and farther north, and only present at latitudes farther south.
We find that the southern spawned eggs from spawning from around late March to mid-April have the highest survival, but also the highest variability between years. This spatial gradient is different from what we would expect as observed spawning is centered around Lofoten (SS 3-5) but indicates that there are additional drivers that shape spawning behavior not addressed in our model. Such drivers could be spawning migration costs (e.g., Jørgensen et al., 2008), spatial patterns in predation pressure (Langangen et al., 2016) and predation avoidance (e.g., Husebø et al., 2009). The seasonal pattern in modeled offspring survival resembles the seasonal pattern in the spawning of NEA cod (peaking around 1st of April, Ottersen et al., 2014), supporting that this timing provides the best balance for offspring to be spawned late enough to benefit from high temperature and fast development but early enough to reach the Barents Sea and FIGURE 8 | Correlation matrix between variables that indicate survival and match-mismatch between the cod larvae 15 days post-hatch within the polygon and their prey of correct size, Calanus finmarchicus (C. fin) nauplii stage III (N3) abundance, larval length and weight at the end of simulation (30th of September) and survival from the egg stage to age 1 (accounting for size-dependent survival through the first winter). Temperature included in this analysis is the mean temperature from March until August within the main distribution area of cod larvae (the polygon in Figure 1) and it is based on the ocean model. DOY corresponds to the day of the year where highest abundance was observed.
have enough food to take advantage of the high growth potential (consistent with, e.g., Kristiansen et al., 2007).

Temperature and Prey, and Their Effects on Growth and Survival
Our findings are consistent with literature indicating that low temperatures result in prolonged stage duration and enhanced mortality independent of prey availability (Pepin et al., 1997;Otterlei et al., 1999). Most larvae are not food-limited within our model simulations. It is mainly the late-and southern-spawned larvae that experience food limitation. The high temperatures these larvae experience imply that the larvae require higher food concentration than other larvae to grow optimally. High temperatures have the potential to secure rapid growth and development (Daewel et al., 2011) and thereby to reduce the number of potential predators, but success requires high food abundance, otherwise high temperatures may be prejudicial (Fouzai et al., 2015). Our results suggest that late-spawned larvae experience higher prey abundance in warm than normal years, but that the food requirements increase more than the food concentrations, leading to increasing food limitations on growth. The food limitations are, however, not so large that the high temperatures are prejudicial for survival.
We find in our model results that the bigger the size at 0age the higher the survival from spawning to age 1. However, the modeled survival is weakly but positively correlated with recruitment success at age 3 and not significantly, so it is not possible to establish a linear relation between modeled survival and recruitment success. Yet, it is known that larger mean size of larvae and juveniles is associated with higher survival through the first winter and recruitment later at age 3 (Ottersen and Loeng, 2000;Stige et al., 2015Stige et al., , 2019. Growing fast and becoming bigger may increase survival because individuals have outgrown their predator diet size range, due to enhanced ability to escape predators, and to the ability to exploit a wider size range of prey (Cushing, 1995). Larger individuals tend to have higher energy reserves and to be more resilient to starvation than smaller individuals, and a larger body size can also be an advantage when dealing with physical stress (Sogard, 1997;Marteinsdottir and Steinarsson, 1998).

Match-Mismatch
We observe a difference in mean time of peak abundance of 15 dph larvae and their suitable prey, but apparently prey availability does not limit survival or growth of most larvae until the 0-age group. In our model we observe that at the date of peak abundance of 15 dph larvae, the maximum concentration of C. finmarchicus is a little below the minimum concentration found by Cushing (1990, page 280), who writes "Above a density of 5 nauplii/L, the cod larvae are well fed." On the other hand, Folkvord (2005) found that most NEA cod larvae in the sea have a size that is near the maximum for their temperature, suggesting that surviving larvae have not suffered food limitation. Our findings are in accordance with Folkvord (2005), and apparently small larvae need less than 5 ind/L to survive without food limitations.
Our results partly support previous studies that have suggested that in warm years, NEA cod larvae grow and survive better than in cold years because of a better temporal match between the first-feeding larvae and their prey (Ellertsen et al., 1987(Ellertsen et al., , 1989Cushing, 1990). Our results are consistent with empirical findings that C. finmarchicus nauplii abundance in spring is higher in warm than normal or cold years (Kvile et al., 2014;Stige et al., 2015), but do not support temperature-driven shifts in the peak timing of the nauplii (Ellertsen et al., 1987(Ellertsen et al., , 1989 and also show large variability in nauplii levels between years with similar temperatures. Our results are also consistent with empirical findings that link high nauplii levels in spring with large size of NEA cod larvae (Ellertsen et al., 1987(Ellertsen et al., , 1989Cushing, 1990;Stige et al., 2015) and between large larval and juvenile size and high survival to later ages (Ottersen and Loeng, 2000;Stige et al., 2015Stige et al., , 2019. Our model results thereby support that the statistical associations between temperature and recruitment through match-mismatch dynamics have a mechanistic basis, but also highlight that these associations are complex and involve temperature-dependent interactions between processes at different spatial and temporal scales.

Model Refinements, Sensitivity to the Results and Future Studies
The current parameterization in the NEA cod IBM includes egg buoyancy but does not differentiate between eggs according to quality and size as a function of e.g., spawning time or spawning adult conditions. An increased egg diameter supports bigger hatched larvae and has positive effects also on yolk-sac larvae (Pepin et al., 1997). However, egg development, time to hatch, size of the larvae at hatching and time of yolk absorption are more significantly affected by temperature than by egg size (Pepin et al., 1997;Steinarsson and Björnsson, 1999). Since biological parameters for the egg and larva development are more affected by temperature than by egg size, we believe our model provides a fair representation of the natural world.
The success of the feeding process and prey selection depend on many parameters such as prey encounter rate, prey, and predator (larvae) size, light, capture success, ingestion rate, prey selection and others. Wind generated turbulence may also affect feeding conditions. Turbulence contributes to the maintenance of the larvae at depths where there is less predation risk than closer to the surface, and that also increases prey encounter rates, suggesting that turbulence can be important to avoid starvation of first-feeding larvae, especially in low prey densities (Sundby, 1997;Kristiansen et al., 2014). On the other hand, turbulence can have a detrimental impact on encounter rates when being too high through the advection of the prey out of the visual range of the larvae (Fiksen and MacKenzie, 2002). Turbulence is not included in the cod larvae IBM used in this study. As food seems to be limiting growth only for the late spawned individuals, we believe that turbulence-enhanced predator-prey contact rates would only have importance for these individuals.
Predation and cannibalism mortality varies in space and time and consequently can affect the spatio-temporal dynamics of survival and future recruitment. These sources of mortality are not included in the models coupled here. The natural mortality rates of NEA cod larvae have been found to vary spatially (Langangen et al., 2014a). In this regard, it is important to understand the spatial heterogeneity of mortality and survival. Age-0 NEA cod has decreased survival when the overlap between juvenile and adult cod increases, likely due to cannibalism (Ciannelli et al., 2007). For the Atlantic cod stock in the North Sea, not only the number of predators is important, but also how they are distributed and at what larval life stage they affect survival rates (Akimova et al., 2019). Akimova et al. (2019) also highlights the importance of the patchiness distribution of both predators and prey.
The predator-prey interaction occurs at much finer spatialtemporal scales in marine ecosystems than the scales modeled here, which can cause uncertainty in predator-prey interactions. Considering the uncertainty in predator-prey and cannibalism interactions mentioned above, it is important to ponder adding further biological features to the early life stages of individuals, as well as further sources of mortality. As topics of future studies, we also suggest analyzing the effect of the spawning adult stock and the variation in demographic characteristics and how this can affect growth and survival of the offspring.

DATA AVAILABILITY STATEMENT
Data can be made available upon reasonable request.

AUTHOR CONTRIBUTIONS
CE contributed with evaluation of model, model analyses, and writing of the manuscript. MS contributed with model runs and model analysis. LS, LC, and FV contributed with constructive discussions, suggestions for model analysis and feedback on the manuscript writing. All authors contributed to the work and contributed to the discussions of hypothesis, methods, and results.