Subsurface maxima in buoyant fish eggs indicate vertical velocity shear and spatially limited spawning grounds

Observed vertical profiles of buoyant particles, in this case pelagic Northeast Arctic (NEA) cod eggs, occasionally deviate from the vertical diffusion‐buoyancy balance by displaying subsurface maxima. Here, we present a mechanism that may explain this phenomenon by combining in situ measurements of NEA cod eggs and concurrent environmental conditions with biophysical modeling of Vestfjorden, Norway. Due to limited observational information, we constructed a spawning season by dispersing eggs with an individual‐based biophysical model forced by a three‐dimensional ocean model including data assimilation improving upper ocean stratification. We show that transient subsurface maxima in eggs are caused by the combination of vertical velocity shear and spatial limitations of spawning grounds. This demonstrates the need for resolving upper ocean small‐scale dynamics in biophysical models to predict horizontal and vertical planktonic dispersal. This is also a precondition for predicting environmental exposure along drift routes, including natural and anthropogenic stressors.

Prediction of particle transport in the upper ocean is needed in a range of applications, from locating accumulation zones of plastic debris (Lebreton et al. 2012), forecasting dispersion of oil spills (e.g., Jones et al. 2016) to predicting distribution patterns of marine planktonic organisms (e.g., Hidalgo et al. 2011). It is useful in survey design for mapping abundances of fish egg, which in turn may be used to estimate the spawning stock biomass as basis for fisheries quotas (Sundby and Bratland 1987;Checkley et al. 1997;Stratoudakis et al. 2006). In case of oil spills, a correct representation of plankton dispersal may be used for quantifications of plankton exposure to oil (Vikebø et al. 2013). Because ocean currents vary with depth, it is necessary to know the dynamical vertical positioning of plankton to obtain correct estimate of dispersal. Hence, simulations of transport and spatial distributions of particles in the ocean must be based on correct representation of physical processes, both horizontal and vertical, from wind-driven currents, fronts and eddies to vertical mixing by wind, tides, and convection. In contrast to dissolved substances that do not influence the specific gravity of the solution (such as inorganic nutrients, dissolved organic carbon, and pollutants of water-based liquids), and hence follow the advection and diffusion of the water masses, particulate matter, such as fish eggs, air bubbles, plastics, and dispersed oil have distinct vertical velocities determined by their buoyancies. Analytical models for the dynamic vertical distribution of particulate matter have been developed for fish eggs (Sundby 1983), air bubbles (Thorpe 1984), and dispersed oil (Paris et al. 2012). Lagrangian biophysical particle-tracking models advect particles according to the current fields of general circulation models (GCMs), while simultaneously adding vertical movements determined by buoyancy (e.g., Ådlandsvik and Sundby 1994;Thygesen and Ådlandsvik 2007) or by vertical behavior in plankton (Vikebø et al. 2007). Variability in ocean currents, horizontal and vertical, caused by atmospheric forcing, tides, or eddies result in particle spreading.
Modeling realistic vertical velocity shear is typically challenging for GCMs, as they tend to smooth out vertical gradients in temperatures and salinity due to inaccuracy in representing turbulence and diapycnal mixing in sigma-coordinate models. This causes, e.g., erroneous flux of energy from the atmosphere and down through the water column. This reduces the accuracy of dispersal forecasts of particles in the ocean, which is particularly relevant near boundaries where stress is exerted (e.g., surface wind stress). For example, Northeast Arctic (NEA) cod eggs and larvae drifting close to the ocean surface are shown to spread more than those deeper down due to the stronger influence of variable meteorological forcing near the surface (Vikebø et al. 2005(Vikebø et al. , 2007.
In this study, we observed vertical distributions of fish eggs and made concurrent measurements of the environmental conditions. We then employed a numerical particle tracking simulation to assess the influence of combined verticalhorizontal processes affecting the vertical distribution of buoyant, drifting particles. In particular, we investigate observations of subsurface maxima in vertical profiles of buoyant NEA cod eggs in the surface mixed layer and show how such a profile structure can be explained.
The choice of NEA cod eggs for this study has several reasons. First, there are sufficiently high concentrations of eggs released in a confined period of time at distinct locations to get excellent observational resolution (Sundby and Bratland 1987). The majority of spawning occurs from March through April (Ellertsen et al. 1989) along the Norwegian coast (Sundby and Nakken 2008). Time of spawning is dependent on the integrated temperature from autumn to spawning time (Kjesbu et al. 2010). The main spawning grounds are in Lofoten, with one spawning hotspot near Henningsvaerstraumen (Sundby and Bratland 1987). Second, more knowledge of the physical processes affecting early life stages of NEA cod is needed to understand the mechanisms regulating survival during the early life stages (Hjort 1914;Ottersen et al. 2014;Strand et al. 2017). Third, the analytical model of vertical distribution of pelagic eggs (Sundby 1983) was based on this species' egg data of physical-biological attributes (i.e., size spectrum and buoyancy) from the present location. Further knowledge about physical-biological attributes throughout egg incubation has been studied in detail by Jung et al. (2014). There are strong vertical and horizontal gradients in egg concentrations at the spawning grounds (Solemdal and Sundby 1981;Sundby and Bratland 1987). The patchy horizontal distribution stems from the spawning behavior (Sundby and Bratland 1987) as well as the presence of physical oceanographic structures. The vertical distribution of NEA cod eggs is determined by the balance between the vertical velocity of the eggs (determined by their buoyancy) and mixed-layer turbulence represented by the vertical eddy diffusivity (e.g., Sundby 1983;Sundby and Kristiansen 2015), see Eq. 1. For pelagic eggs (i.e., eggs with density lower than the density of the upper mixed layer), such as NEA eggs, concentration declines exponentially with depth from the surface in proportionality to egg ascending speed and in inverse proportion to the eddy diffusion coefficient (Sundby 1983(Sundby , 1991. Still, observed vertical profiles of NEA cod eggs occasionally reveal a transient subsurface maximum inconsistent with the steady-state vertical analytical formulations by Sundby (1983) (see, e.g., Sundby 1983;Röhrs et al. 2014). Here, we investigate potential mechanisms causing these observed deviations from the analytical vertical model for pelagic fish eggs.
Our main hypothesis is that subsurface maxima in buoyant particles may occur due to horizontal movement in the presence of vertical velocity shear and strong gradients in horizontal egg concentration, conditions that are often observed at spawning hot spots. By designing a numerical experiment based on measured conditions at the main spawning ground of NEA cod, we quantify the frequencies by which such events occur, and explore which conditions favor such incidents. We then look at how adding data assimilation in the GCM can improve the representation of vertical and horizontal shear and compare with observations of vertical egg distribution.

Materials and methods
In our analysis of the main hypothesis, we initiate particles at a well-known spawning site inside Vestfjorden and model their subsequent dispersal with an individual-based biophysical particletracking model (Ådlandsvik and Sundby 1994) forced both by idealized currents and hourly current velocity, hydrography and turbulence from a three-dimensional GCM, constructed by the use of state-of-the-art data assimilation methods (Sperrevik et al. 2017). Particles are initiated at multiple adjacent point locations so that we may analyze the effect of narrow vs. wide spawning grounds for the occurrences of subsurface maxima. An evaluation of the modeled ocean circulation is given in Supporting Information Fig. S1 where we compare two model realizations, with and without data assimilation, against in situ measurements of NEA cod eggs from observations in 1984 (Sundby and Bratland 1987). By including data assimilation, the NEA cod dispersal is improved compared to the observations (Supporting Information Fig. S1). The year 1984 is chosen because of the extensive measurement campaigns providing both physical data for assimilation and evaluation of the GCM, and egg distribution data. As only horizontal (and not vertical) egg observations were available in 1984, we compare our modeled vertical egg profiles with corresponding observations from a scientific cruise during the spawning season from 4-7 th April 2016 in the same area. Then, only vertical egg profiles were sampled, but no horizontal egg coverage were carried out. However, egg data from 1984 and 2016 may be compared because the salinity in the area has not changed enough to affect the physical NEA cod egg buoyancy differently. See "NEA cod eggs as oceanographic drifters" section below. The horizontal patchiness of the eggs observed in 1984 mirrors the typical traditional spawning areas repeatedly observed during earlier studies (Ellertsen et al. 1981b(Ellertsen et al. , 1984Sundby and Fossum 1990;Sundby et al. 1994), including the specific area sampled in 2016, the Henningsvaerstraumen. This was studied in detail by Sundby and Fossum (1990), where typically horizontal egg concentration decreases by two orders of magnitude over a 10 km distance from the center of the spawning area (Sundby and Bratland 1987).

NEA cod eggs as oceanographic drifters
An important characteristic of marine fish eggs is that they are homohaline, implying that they by osmoregulation maintain Strand et al. Subsurface maxima in buoyant fish eggs constant internal salinity independent of the ambient salinity (Sundnes et al. 1965). Furthermore, fish eggs are ectotherms, meaning that their internal temperature equals that of the environment. The thermal expansion coefficient of the eggs is approximately equal to that of the ambient seawater (Sundby and Kristiansen 2015), allowing the in situ buoyancy to be calculated through laboratory experiments based on salinity alone and independent of the ambient temperature. Hence, the buoyancy of NEA cod eggs depends on salinity, but not the temperature. The laboratory-based neutral buoyancy of NEA cod eggs expressed in salinity units ranges between 29.5 psu and 33.0 psu, with an average neutral buoyancy of about 31.0 psu (Solemdal and Sundby 1981 S2). This means that a small portion of the observed eggs, in both years, could be denser than the ambient water masses. Based on data from Jung et al. (2012a,b) and Stenevik et al. (2008), this amounts to 5.1 (4.1)% of the eggs for 32.7 (32.8) psu, giving a difference of only 1% in potentially denser eggs between 1984 and 2016. A second important characteristic is that the NEA cod at the spawning areas in Lofoten release their eggs in the thermocline, within a temperature range of 4-6 C, usually at depths varying between 50 m and 200 m (Ellertsen et al. 1981a). The thermocline defines the interface between the Norwegian coastal waters (cold and relatively fresh) and inflowing Atlantic waters (warmer and more saline), a typical hydrographic situation for the Lofoten spawning area during spring time (Ellertsen et al. 1981b). While the upper ocean temperature is higher in 2016 than in 1984, the spawning still occurs in the transition zone between Atlantic and Coastal waters.
These characteristics make NEA cod eggs positively buoyant, and newly spawned eggs rise toward the surface and reach their equilibrium vertical distribution in less than 24 h. The exact time to equilibrium depends on the intensity of wind mixing (Sundby 1991). At the ambient temperatures of upper layers of the Lofoten spawning areas, NEA cod eggs typically hatch after about 3 weeks (Strømme 1977), allowing considerable drift distances in the upper ocean from the spawning grounds toward the nursery area before hatching into the larval stages. NEA cod egg develop through six defined stages (Strømme 1977), enabling quantification of how long individual eggs found at sea have been adrift.

Observations from a scientific cruise, 4-7 th April 2016
A scientific cruise was conducted 4-7 th April 2016 with R/V Johan Hjort by the Norwegian Institute of Marine Research in collaboration with the Norwegian Meteorological Institute, the Nansen Environmental and Remote Sensing Center, and the National University of Ireland, Galway. Nine vertical egg profiles (see Fig. 1 for locations), including egg-stage determination according to Strømme (1977), are used together with the oceanographic and meteorological observations to evaluate and compare with our modeling study, as explained below.
Eggs were sampled with a Xylem submersible electric pump with a pump capacity of about 100 L min −1 . Seawater was pumped on deck through a 75-mm hose and filtered through a T-80 plankton net with mesh size 375 μm. Pump samples were taken at 1 m depth (except for one location at 1.5 m), then every 5 m from 5 m to 30 m. Sampling volume from each depth was 200 L. Egg profiles are presented in numbers m −3 . The measurement increment is 1 egg per 200 L, i.e., 5 eggs m −3 . The pump technique of sampling vertical egg profiles is therefore sensitive when low numbers of eggs are observed. Given the actual sea state during the cruise, there is an assumed uncertainty of 0.25 m per measurement depth due to the movement of the ship. In addition to the vertical egg samples, 39 net hauls were sampled with the same plankton net. In total, from both vertical egg samples and net hauls, 2991 eggs were counted and staged. For every egg profile, a CTD profile was taken with a Seabird CTD instrument (SBE 911+). The data was postcalibrated against water samples taken with every CTD profile.
Current velocities were measured using an acoustic upward looking Aanderaa Recording Current Doppler Profiler (600 kHz) placed at 40 m depth on a mooring in the center of the survey area (68.09 N, 14.07 E) which is at the assumed center of the Henningsvaerstraumen spawning area. The bottom depth at the mooring location is 105 m. The instrument was operational from 4 th April 13:11 UTC until 7 th April 03:46 UTC. The upper 5 m of the data before 5 th April 10:00 UTC could not be used due to higher frequency interference with another instrument working at the beginning of the cruise. Processed velocity data are stored in 5-min averages in 1 m depth intervals. The measurements are filtered with a Hanning window to remove variability of time scales less than 1 h to obtain the same time step as the other observations. Automatic wind observations were taken from the nearest meteorological station Skrova Lighthouse (WMO st.no. 01160), located on a small island 11 m above sea level (68.15 N, 14.65 E) 25 km from the observational site and operated by the Norwegian Meteorological Institute (http:// eklima.met.no, accessed 6 th April 2017). A comparison with the wind mast from the ship shows similar observations, though a time lag of a few hours depending on the situation and location of the ship.

Ocean model setup
Particles are transported by hourly three-dimensional current fields from a Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams 2005;Haidvogel et al. 2008). We used the Generic Length Scale mixing scheme with k-ω setup for  Beldring et al. 2003), and eight tidal constituents from the TPXO global inverse barotropic model (Egbert and Erofeeva 2002). Initial and boundary conditions are from the SVIM hindcast archive with a 4 km × 4 km horizontal resolution (Lien et al. 2014). A reanalysis of the ocean circulation was produced by the use of four-dimensional variational (4D-Var) data assimilation (Sperrevik et al. 2017) using hydrographical observations from an extensive field campaign performed by IMR in the main spawning area for NEA cod, the Vestfjorden, in 1984 (Sundby and Bratland 1987) as well as satellite sea surface temperature.
Individual-based particle tracking NEA cod eggs are released continuously at point locations in a regular grid, centered around the main spawning ground at Henningsvaerstraumen (Sundby and Bratland 1987). The eggs are advected hourly using a 4 th order Runge-Kutta advection scheme. The model variables are tri-linearly interpolated to the individual time-varying locations of each egg. The buoyancies of eggs are based on the individual egg sizes and densities (see Eq. 1 and Sundby 1983) and modeled ocean densities. Vertical dynamical positioning of eggs is calculated based on the numerical scheme by Thygesen and Ådlandsvik (2007) utilizing the turbulence from ROMS at the individual time-varying location of each egg (see, e.g., Röhrs et al. 2014). The spawning ground is represented by 66 locations (one location per grid cell inside the box, see Fig. 1). Particles are released at 50 m depth, with 25 particles per location every 6 h for 60 d, corresponding to the main spawning period from 1 st March to 30 th April, resulting in a total release of 397,650 eggs. The individual-based biophysical particle-tracking model is run for 80 d, ending 20 th of May to ensure that all initialized eggs have hatched. The eggs mature and hatch according to ambient water temperature (Folkvord 2007). Particle positions are stored every 3 h, resulting in eight track positions per day, thus resolving tidal motion.

Sensitivity analysis of subsurface maxima under idealized currents
A sensitivity analysis is included where eggs are transported by currents resulting from a two-step reduction of the original modeled currents to an artificial constant depth-independent horizontal current equal to 1 cm s −1 . Also, we have tested the importance of dynamical vertical positioning of eggs due to turbulence and buoyancy for subsurface maxima by adding simulations with constant egg rise velocities of 1 mm s −1 .

Sampling vertical profiles of cod eggs in the model simulation
Cod eggs in Lofoten hatch after about 3 weeks. Larvae have different buoyancy than eggs as well as having a vertical behavior. Therefore, the particles are removed from the model once they hatch. Vertical profiles of eggs in the model simulation are then sampled to search for subsurface maxima, in the upper 20 m where the majority of eggs are located. To identify subsurface maxima, vertical profiles of particles are sampled at every grid cell where particles were initiated (in the spawning ground set to 66 grid cells). Two methods of sampling vertical profiles of particle distributions were carried out: (1) by only considering the particles initiated at the same grid cell as they are subsequently being sampled and (2) by considering all particles independently of where they are initiated. The latter is expected to better reflect the observations with the egg-pump stations in 2016 as spawning cod is not all gathered in a single point location. The comparison of the two ways of sampling the model enables us to consider the effect of vertical shear and spatial extent of the spawning ground on vertical profiles of buoyant eggs in a turbulent environment, particularly whether situations with subsurface maxima initiated by vertical shear are obscured by horizontal transport between neighboring spawning locations. Since the biophysical model is computationally demanding, there are a limited number of particles representing eggs. We therefore test the robustness of our results by comparing with occurrences of subsurface maxima requiring a minimum amount of eggs present. The minimum threshold is tested for values of 30%, 40%, and 50% of the mean surface egg concentration.

Analytical vertical profile of NEA cod eggs
The mean vertical egg concentration C(z), from a balance between turbulent mixing and buoyancy of the NEA cod eggs, is given from Eq. 5 by Sundby (1983): where z is depth, C a is the known egg concentration at depth a, K is the eddy diffusion coefficient, and w is the ascending velocity of the eggs.
Here, the following values for the variables have been used to calculate a mean vertical egg profile: Depth a = 1 m (the surface layer of the model), K = 0.02 m 2 s −1 according to Eq. 18 in Sundby (1983) with mean wind speed of 7.4 m s −1 calculated from NORA10 March-May wind, w = 1 mm s −1 from Fig. 1 in Sundby (1983) with mean NEA cod egg diameter of 1.4 mm (from Solemdal and Sundby 1981) and the density difference (Δρ) between the ambient water and NEA cod eggs is 1.8 kg m −3 (ρ water − ρ egg = 1026.6-1024.8 kg m −3 ).

Observations 4-7 th April 2016
Variations in vertical NEA cod egg profiles, including subsurface maxima, are observed at Henningsvaerstraumen, one of the main spawning areas of NEA cod during the cruise 4-7 th April 2016 (Fig. 2). The concurrent mean water salinity increases almost linearly with depth from 33.0 (range: 32.7-33.3) psu at the surface to 33.3 (range: 33.1-33.5) psu at 30 m (Supporting Information Fig. S2). Subsurface maxima (Fig. 2a) occur during periods of enhanced northeasterly wind (Fig. 2b,c).
The mean current during the cruise is 0.15 m s −1 (Fig. 2c) which corresponds to a displacement of about 65 km in 5 d. Observed ocean currents at multiple depths display vertical current shear (Fig. 2d,e). During the calm wind period succeeding a strong south-westerly wind event (4 th April 13:00 UTC to 5 th April 19:00 UTC, Fig. 2), the ocean current speeds are below 0.15 m s −1 and generally increasing with depth bearing southeast. Subsequently, the wind strengthens and veers north-easterly (5 th April 19:00 UTC to 6 th April 16:00 UTC, Fig. 2) with stronger ocean currents to the south, particularly near the surface. Finally, the winds weaken while maintaining bearing, though the ocean currents turn northeasterly (6 th April 16:00 UTC to 7 th April 08:00 UTC, Fig. 2).

Numerical model 1 st March 1984-20 th May 1984
As particles are being released, they rapidly adjust to the modeled ambient density structure and vertical mixing, resulting in profiles with near exponential decrease from the surface to about 30 m (Fig. 3). Figure 3a shows all vertical profiles sampled by method 2 (considering all particles independently of where they are initiated), every 3 h, at grid cell 28 (approximately the center grid cell, Fig. 1), revealing that while the median profile decreases from just below 100 eggs m −3 at the surface to almost none at 30 m depth, there are incidents when the surface concentrations are an order of magnitude higher. Figure 3b shows that vertical profiles vary significantly depending on where, within the modeled spawning ground, they are sampled. Profiles at grid cell 28 of the modeled spawning ground on 17 th April are distinctively different from those at the southern or northern boundary about 12 km away. The time evolution of vertical egg concentration at grid cell 28 shows a large vertical variability with distinct periods of enhanced mixing and eggs distributed deeper, e.g., early in March (Supporting Information Fig. S3). The focus in this study, however, is on the occurrence of subsurface maxima (Fig. 3c). Figure 4a resolves occurrences of subsurface maxima in the upper 20 m of the vertical egg profiles through time and space. In general, the southern grid cells have the most frequent occurrences of subsurface maxima in egg concentrations (gray dots), while the northern grid cells have the most frequent occurrences if adding the threshold of 40% described in "Materials and methods" (red dots). Summarized across all grid cells per time step (Fig. 4b), the time series show periods with increased occurrences of subsurface maxima, and six shorter periods where none of the grid cells have subsurface maxima (9 th and 21 st March, 6-7 th April, 11-13 th April, 1 st -3 rd May, and 18 th May). Only counting the particles spawned locally (method 1), thereby not taking into account import of particles from neighboring grid cells, the occurrence of subsurface maxima increases substantially (Supporting Information Fig. S5, without including the threshold). On average over the whole period (1 st March to 20 th May), 38% of the area of the spawning ground has subsurface maxima when not allowing import of particles (Supporting Information Fig. S5), while 22% of the spawning ground has subsurface maxima when allowing import of particles (gray dots in Fig. 4a,b). The latter number reduces to 10/8/6% if adding the egg threshold requirement of 30/40/50% while all have similar spatial variability (red dots, Fig. 4a,b). In order to investigate the causes of this variability, concurrent time series of currents at three depths (the spatial mean across the 66 grid cells) and wind at grid cell 28 are analyzed (Fig. 4c,d).
The sum of occurrences of subsurface maxima in vertical profiles of egg concentrations (without including the threshold) is correlated with the wind speed and the surface current, where the Pearson linear correlation coefficient r is 0.63 (p < 0.001, t = 20.3, df = 639) for wind and 0.46 for surface current (p < 0.001, t = 13.2, df = 639). Correlating wind speed directly with the sea surface current gives r = 0.73 (p < 0.001, t = 27.4, df = 639). Testing for time lags between subsurface maxima and forcing do not result in significant improvements of the correlations (wind: r = 0.65, surface current: r = 0.47). Replacing surface current with current shear represented as the difference between 20 m and surface, or 7 m and surface, results in about the same  Subsurface maxima do not necessarily occur simultaneously throughout all 66 spawning grounds (Fig. 4a). Focusing on 20 th -30 th March, there are at first no subsurface maxima (21 st March) then subsurface maxima throughout the spawning ground (gray dots, Fig. 4a), with a northward subsurface maximum signal propagating through the spawning ground (red dots, Fig. 4a). A progressive vector diagram for grid cell 28 of the same period displays a strong current shear with a varying direction and decreasing strength with depth ( Fig. 5a). Red dots (Fig. 4a) show that a collection of many eggs is advected across the spawning ground but that the near surface ones are continuously shed off due to the shear resulting in subsurface maxima. Focusing on a second period, 15 th -30 th April, there are two distinct periods of enhanced currents above 0.25 m s −1 , corresponding to similar peaks in wind forcing and subsurface maxima. Also, the corresponding progressive vector diagram at grid cell 28 shows a vertical velocity shear of decreasing current strength with increasing depths (Fig. 5b). Apparent, in Fig. 5b, there is also a two-layer stratification with coastal waters on top (here; upper~25 m), and Atlantic water below (also seen in lower left panel, Supporting Information Fig. S1). Subsurface maxima appear first along the southern rim of the spawning ground and later also to the north. This suggest that despite that the velocity shear may result in subsurface maxima, this is delayed to the north because of continuous supply of eggs near the surface from upstream sources. As the source empties, subsurface maxima appear successively northward.
Reducing the current shear or introducing a fixed rise velocity of eggs reduces the number of profiles with subsurface maxima at grid cell 28 (Supporting Information Fig. S4, only considering eggs initialized at the same grid cell, method 1). The left panel shows the original number of total profiles between 1 st March and 20 th May with subsurface maxima (262 profiles). The mid panel shows how this number decreases (237 profiles) if moving particles with reduced current velocity shear. Right panel shows further reduction to about half the number of occurrences if only using the constant current (129 profiles). Removing turbulent dynamical vertical positioning of eggs (not shown) and combining with either the modeled currents or the constant current result in either a strong reduction or a complete removal of subsurface maxima.

Discussion
Dispersal of buoyant particles depends on their vertical positions and the vertical current shear. Theoretical considerations of vertical distribution provide a mathematical framework for quantifying vertical profiles under various oceanographic conditions, given their individual densities and sizes of the particles (e.g., Sundby 1983;Thorpe 1984). Concentrations of such particles decrease exponentially with depth, where the vertical gradient depends on the buoyancy of the particles and the ambient level of turbulence in the water column. Occasionally, we measure vertical profiles of buoyant particles in the field, in this case NEA cod eggs, that differ from the vertical diffusion-buoyancy balance and instead display subsurface maxima. This is not because theory is proven wrong, but because additional horizontal processes are interacting.
Our main hypothesis is that the deviations in the vertical profiles from the diffusion-buoyancy balance are caused by the combination of vertical velocity shear and strong horizontal gradients in egg concentrations around a spawning ground. Both conditions are observed by extensive measurement campaigns, as reported here and previously (Sundby and Bratland 1987). Strong horizontal gradients and the presence of older eggs therefore support previous findings that Henningsvaerstraumen is characteristically a retention area compared to spawning grounds outside Lofoten (Sundby and Bratland 1987), though older eggs may also originate from spawning grounds elsewhere. Since there is sparse spatial information from observational cruises to analyze, we constructed a spawning season where we perform an extensive egg survey in a numerical model. We demonstrated that there is large variability in the egg profiles, with subsurface maxima occurring transiently during periods of otherwise exponential decaying concentrations of eggs.
The mean occurrence of subsurface maxima is 22% over the whole period when accounting for import of eggs from neighboring spawning grounds. If only accounting for eggs spawned at the site sampled, the mean occurrence increases to 38% illustrating that spawning grounds with a limited horizontal extent have a higher propensity of exhibiting subsurface maxima in egg concentration. Including the egg concentration threshold illustrates that a signal of subsurface maxima should be treated with care if there are low numbers of egg sampled. By investigating periods of increased occurrences of subsurface maxima against the wind forcing and ocean currents, we find that both factors favor subsurface maxima, in particular periods of persistent forcing or after sudden transitions in the direction. We find a significant correlation between the time series of wind speed or surface current speed and the occurrence of subsurface maxima. Considering time lags between wind, currents and subsurface maxima did not give significant improvement but indicates that wind and currents may lead subsurface maxima by a few hours. A model study with higher temporal resolution would enable a more decisive answer. There is lower correlation with subsurface maxima against surface current (r = 0.46) than subsurface maxima and wind speed (r = 0.63), which seems counter-intuitive considering that the ocean current is the direct forcing. However, the wind represents the transient energy exerted on a relatively large area affecting vertical distribution of eggs through several physical processes including ocean currents and turbulence. Contrary, measured and modeled ocean currents represent dispersion and shear on a scale smaller than the size of the spawning ground here represented by 66 grid cells. Hence, even within these spawning cells, the currents vary, and the mean current across the spawning ground therefore correlate less with subsurface maxima than with wind.
A sensitivity test shows that with a constant horizontal flow and no vertical dynamical positioning due to turbulence, there will be no subsurface maxima. However, increasing the strength of the horizontal flow or reducing the width of the sampled water column will eventually result in subsurface maxima because eggs are moved outside the sampling area before they surface. In reality, turbulence opposes this as it contributes to erase vertical gradients introduced by spawning at depth and join forces with buoyancy moving eggs toward the surface.

Other processes potentially causing subsurface maxima
A vertically varying eddy diffusivity coefficient K(z) may affect the rate at which particles are redistributed vertically if introduced at a certain depth but cannot cause a vertical gradient of particle concentrations unless the particle density is lower than the ambient water density (e.g., Thygesen and Ådlandsvik 2007). Hence, a high level of turbulence near the sea surface cannot cause egg aggregation immediately below but combined with low numbers of eggs, this may happen by chance because of their turbulence-induced dynamical induced distribution. This is supported by the sensitivity analyses quantifying occurrences of subsurface maxima with and without turbulence. Increased wind forcing results in increased vertical mixing in the ocean causing buoyant particles to be mixed down through the water column. Compared to calm conditions, the concentrations of eggs still decrease near exponentially with depth, but to a much lower degree. As the wind forcing dies off and the mixing level ceases, buoyant eggs start to rise toward the surface. For example, NEA cod eggs in Henningsvaerstraumen mixed down to 10 m depth will rise toward the surface within the next $ 3 h, assuming a typical ascending velocity of 1 mm s −1 . If a strong salinity structure re-establishes before the eggs have reached a new vertical profile, eggs rise faster at depth than near surface because the buoyancy decreases. Altogether, this may cause transient subsurface maxima in eggs.
Egg densities vary, as mentioned in "NEA cod eggs as oceanographic drifters" section, where the neutral buoyancy expressed in salinities ranges from 29.5 psu to 33.0 psu (Solemdal and Sundby 1981). From CTD profiles in 2016 (Supporting Information Fig. S2), the mean salinity observed in the surface layer is between 32.7 psu and 33.3 psu, which is the upper neutral buoyancy range of the NEA cod eggs. This makes 5.1% of the eggs potentially negatively buoyant (heavier than the ambient water masses) and able to sink creating subsurface maxima. The two lightest salinity profiles correspond to the egg profiles #304 and #305. Since there is low number of eggs sampled here, these measurements may be sensitive to the pump measurement technique. At profile #306, #310, and #315, however, there are higher numbers of eggs sampled experiencing subsurface maxima. Since the salinity is above 33.1 psu for these profiles, the subsurface maxima signals cannot be explained by buoyancy differences in the eggs, and other physical processes must be responsible.
Newly spawned eggs will attain a vertical equilibrium distribution within about 24 h, depending on the induced mixing (using Eq. 1, see Sundby 1983;Sundby 1991). Hence, a subsurface maximum could be observed for newly spawned eggs (Sundby 1991), but only for a short time period in the beginning of the spawning season making it an unlikely process to be observed.

Strand et al.
Subsurface maxima in buoyant fish eggs Langmuir cells created by Stokes drift (induced by the presence of surface waves) potentially cause inhomogeneous mixing (Grant and Belcher 2009;Belcher et al. 2012;Harcourt 2013), but only in narrow bands of the ocean making it unlikely that they will actually be observed during the time period it takes to measure the vertical egg profiles.
Air bubbles are introduced into the upper ocean when waves break, giving a theoretical possibility, if the air bubbles are small enough and high enough in numbers, to affect the density of the water column. This effect is confined to the upper few meters, on a vertical scale comparable to the significant wave height (Scanlon et al. 2016). However, upper concentration estimations of air bubbles with diameter > 10 μm are 10 6 m −3 (Zhang et al. 2010). With cod eggs of 1 mm, this gives approximately 1/1000 bubble per egg volume which is too low to change the buoyancy of the individual eggs.
Selective predation on NEA cod eggs in Vestfjorden could be by planktonic predators in the upper water column, with Norwegian Spring-Spawning herring and jellyfish being potential candidates during this time of the year (March-April). The massive numbers of predators needed in the upper 5 m of the water column in order to reduce the egg concentration significantly in a short enough period of time makes this unlikely.
The available turbulence schemes in ROMS are forced by the boundary condition to give turbulent diffusivities that approach zero at the surface Röhrs et al. 2014). This gives a bias in the upper layer compared to the analytical solution as diffusivity is underestimated. In turn, this causes too many cod eggs at or near the surface counteracting the presence of subsurface maxima in the model as compared to observations. Sperrevik et al. (2017) found that with data assimilation in the GCM, the water column becomes more stratified so that a shallower part of the water column responds to wind forcing exerted at the surface. Again, a higher spatial resolution in the GCM both vertically and horizontally would provide added details on the dynamical manifestations of wind forcing on the upper ocean.

Conclusion
Observations from a cruise in Lofoten, Norway, in 2016 reveal transient subsurface maxima in NEA cod egg profiles consistent with previous observations. Our main hypothesis is that this is caused by spatially limited spawning grounds and the presence of vertical current shear. By running a high-resolution model with assimilation of available hydrographic data, we were able to reproduce subsurface maxima and relate this to wind stress and vertical current shear. An idealized sensitivity analysis shows that if vertical current shear is gradually reduced for egg dispersal from a spatially limited spawning ground, then the occurrence of subsurface maxima also decays and eventually disappear, supporting our main hypothesis.