Risk of oil contamination of fish eggs and larvae under different oceanic and weather conditions


 An oil drift model is applied to determine the spread of oil spills from different locations along ship lanes off southern Norway every month for 20 years. These results are combined with results from an egg- and larvae drift model for Atlantic cod (Gadus morhua) to determine their risk of being impacted by oil. The number of eggs and larvae exposed to oil contamination is connected to environmental conditions. The highest risk of overlap between an oil spill and cod in early life stages occurs during March and April when the eggs and larvae concentrations are highest. Spills off the west coast pose a greater risk because of the ship lanes’ proximity to the spawning grounds, but there is large interannual variability. For some spill locations the interannual variability can be explained by variability in wind and ocean currents. Simultaneously occurring onshore transports lead to a high-risk situation because both oil and larvae are concentrated towards the coast. This study demonstrates how results from oil drift and biological models can be combined to estimate the risks of oil contamination for marine organisms, based on the location and timing of the oil spill, weather/ocean conditions, and knowledge of the organisms’ life cycle.


Introduction
Marine oil spills can affect the ecosystem both through their physical influence on the organisms and through the effect of the various chemicals in the oil. On the physical side, organisms can get caught in the viscous oil. Oil can also form a barrier to the surface, which can be damaging for some species of fish when they come to the surface to fill the swim bladder (Sundby et al., 2013). There are also several chemicals in the oil that can increase the mortality for eggs and larvae (Neff et al., 2000;Barron et al., 2004). Oil can stimulate the formation of marine snow (Passow et al., 2012) that may subsequently sinks to the ocean floor where it can impact the benthic fauna. If conditions become anoxic, oil in the sediment would decay more slowly since the biological remineralization of the oil is retarded in anoxic conditions. The use of dispersants to increase the rate of remineralization has also received attention in terms of its potential negative effect on the ecosystem (Kleindienst et al., 2015;Vikebø et al., 2015).
This study focuses on the seas around southern Norway (Figure 1), where the currents are dominated by the Norwegian Coastal Current (NCC) that flows along the coast all the way from the Skagerrak to northern Norway. This current would primarily transport oil spills, eggs, and larvae along the coast. In the northern part of the area, the NCC flows side by side with the Norwegian Atlantic Current (NAC), which is a more saline current also directed northwards. Strong winds of variable directions and high waves occur frequently, particularly during winter, and contribute to the upper layer drift. The area has strong mesoscale activity, particularly in the frontal region between the fresher NCC and the eastern branch of the NAC, and thus it is difficult to predict the local currents on short time scales. The tidal V C International Council for the Exploration of the Sea 2019. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. amplitudes, and consequently the tidal currents, are small (<1.0 m) in the southern parts of Norway, but become larger in northern Norway.
The Norwegian Sea and the North Sea are busy shipping areas, in connection with both general cargo transport and oil production. In the last two decades the Norwegian Coast has experienced several large oil spill incidents involving ships, for example the Server-accident in 2007 that contaminated large parts of the west coast (Boitsov et al., 2012).
The currents, winds, and waves vary both interannually and seasonally. Since all these factors will influence the trajectories taken by both an oil spill, fish eggs and larvae it is likely that the risk of an oil spill impacting fish larvae will vary in the same way. The main objective of this study is to determine the seasonal and interannual variability of how much oil spills at selected locations can impact cod eggs and larvae. A further aim is to assess the environmental forcing factors that influences this variability. Although the selections of the spill locations are guided by ship traffic density maps, the spills we operate with are still hypothetical. The actual mortality effect of the oil is also highly uncertain; therefore, we can only estimate the amount of overlap between areas with oil and areas with eggs and larvae, but not quantify the effect in terms of mortality on the eggs and larvae.
Beside the environmental factors, the properties of the oil and the eggs and larvae themselves will also influence the results. The chemical and physical properties of the oil will determine, among others, how easily it can be mixed down in the water column, how fast it will evaporate, how much it will spread and to what degree it will form an emulsion (Reed et al., 1999). Fish eggs can have different positions in the water columns, some species having eggs that sink or stick to the bottom, while others drift at intermediate depths. Cod have positively buoyant eggs and they will therefore be located close to the surface and are vulnerable to oil spills. As they enter the larval stage and the larvae grow, they develop the ability to position themselves vertically in the water column. The vertical position will determine whether they come in contact with a surface oil spill and, because currents can change with depth, how they will drift relative to the oil spill (Vikebø et al., 2007).
While atmospheric reanalysis products have existed for many years (Dee et al., 2011), oceanic reanalysis products have only recently been developed (Sakov et al., 2012;Balmaseda et al., 2015). By assimilating observations of sea-surface temperature and height, sea ice, as well as vertical profiles of temperature and salinity (e.g. from ARGO), ocean reanalyses provide more realistic estimates of the physical model fields compared to free (non-assimilating) model simulations (Lien et al., 2016;Xie et al., 2017). In this study, we use both oceanic and atmospheric reanalysis products to force drift models for oil, eggs, and larvae. The results are used to compute potential impacts of an oil spill at different locations on cod eggs and larvae both the climatological seasonal variations and in the individual years of the simulation period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). While the subject of impact of oil spills on fish larvae by combining oil spill and larvae drift models has been the focus of several previous studies (Spaulding et al., 1983;Vikebø et al., 2014), the novelty in this study is that we used a longer simulation allowing for the evaluation of both climatological and interannual variability of the risk based on analysis of the 20-years of simulation. Analysis of the model forcing fields also makes it possible to determine the environmental conditions favourable for high-or low-risk situations.

The oil drift model
The oil drift model, OD3D, (Martinsen, 1994;Wettre et al., 2001) is a Lagrangian model that takes into account three-dimensional currents to handle oil spills released at depth as well as oil that has been mixed down in the water column. The Lagrangian numerical particles represent parcels of oil (Delvigne and Sweeney, 1988), and the properties of each oil parcel are tracked through time and changes depending on how long it has been in the water (age) and the environmental conditions it is expose to. OD3D has a database containing around 60 types of oils, ranging from light to heavy oils and from different oil and gas fields in the North Sea. These oil types differ from each other with respect to chemical characteristics and how they are affected by various weathering processes such as evaporation, emulsification, and vertical dispersion. For these model runs we have chosen to use "heavy oil" in the spill simulations because it is the most used fuel in commercial shipping in the simulation period, although lighter fuels are becoming more used. Typical for a heavy oil is that there is little or no evaporation or dissolution, it weathers slowly and can persist as floating or beached oil for a week or more. Since it persists longer in the water, a heavy oil can have severe impacts on the marine ecosystem.
The oil drift model requires information about bathymetry and currents, temperature, and salinity at multiple depths; this information is obtained from an ocean reanalysis (see Table 1 and the "Ocean reanalysis and derived fields" section for more details). The oil drift model also needs wind velocity at 10 m ( period (Table 3). Since the forcing comes from different sources, both the spatial and temporal resolutions of each forcing differ (Tables 1-3). The temporal resolution of the geophysical forcing data is daily for the oceanic variables, 6-hourly for the atmospheric variables, and 3-hourly for the wave-related variables. The oil is released as numerical particles, each representing a certain fraction of the total volume of oil spilled. In the remainder of the paper the term "oil particle" refers to these numerical particles.
To assess the oil drift in a certain month, oil is released for 2 h every day over the course of a month at each predetermined location (see "Spill locations" section). This is done to yield a good temporal spread of the release throughout the month, but limit the number of numerical particles to a manageable size. The oil particles are released between midnight and 2:00 am each day, but the time of day of the release should not impact the results since tides are not included in the forcing. At each spill location, 20 particles are released every day, in total $600 over the course of a month. Sensitivity runs using twice the number of numerical oil particles have been performed at a few of the spill locations and there were only small differences in the resulting risk maps described in the "The statistical risk calculation" section. We therefore decided to continue with 20 particles per day per oil spill site. The oil drift model differentiates between the six states of the oil, which are "at the surface," "dispersed," "aged," "evaporated," "permanently submerged," and "stranded" oil. Each oil particle has a maximum residence time of 6 days, after which the particle is considered as aged and removed from the simulation. Particles are no longer tracked by the simulation when they beach, evaporate or are permanently submerged in the water column. Dispersed oil is oil droplets that are below the surface at some depth. These droplets are tracked by the simulation and may resurface under the right environmental conditions to form new slicks or rejoin the original slick. The oil drift model's internal time step is 15 min.

The egg and larvae drift model
We chose cod as our target species. Cod is one of the most valuable commercial fishes in Norwegian waters and, since the eggs are buoyant they are vulnerable to oil spills. The larvae drift model is structured as an individual-based model (IBM: Grimm and Railsback, 2005) that is coupled to the physical environment using a Lagrangian transport model. The particle drift is computed offline, meaning that the particle drift is computed based on stored ocean circulation model fields. The particle transport model is based on a version of ECOSMO-IBM previously used for brown shrimp in the North Sea  and was modified for the use of Atlantic cod in Norwegian waters. The main modifications are that the temperaturedependent stage duration function is adjusted to represent cod development and the buoyancy formulation for the eggs is modified to account for cod eggs being approximately spherical, in contrast to shrimp, which are more elongated. The variables required to force the IBM are three-dimensional current velocities, temperature, salinity, and vertical diffusivity coefficients, all of which are either taken directly from the previously mentioned ocean reanalysis (Table 1) or derived as described in the "Ocean reanalysis and derived fields" section. The internal time-step of the IBM is 30 min.
As in the original IBM, the sum of advective velocity V a and diffusive velocity V d determines the particle displacement at a given time and location: The input from the physical model is linearly interpolated to the particle locations. The vertical diffusion is assumed to be a random process. With the diffusion coefficient (D), calculated by the hydrodynamic model, the diffusive velocity is derived using a Monte-Carlo method. Based on an approach proposed by Bork and Maier-Reimer (1978), we formulated the diffusive velocity as a simple random walk so that V d t; X ð Þ ¼ RÃ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 6D=Dt p , where R½À1; 1 is a random number. In the frame of our experiment cod eggs cannot be considered to be passive particles so that, in addition to the advective and diffusive velocities, they are given a velocity relative to their environment.
As Atlantic cod eggs are known to have a specific density (Jung et al., 2012) an additional buoyancy term must be added to the vertical particle displacement. As in  the formulation is based on solving the vertical component of the equation of motion, where gravity, Stokes friction, and the buoyancy term that results from the density difference between eggs and their environment (Archimedes' principle) are considered. This term is expressed as Equation (2), and discretized in time [Equation (3)]: here w is the vertical velocity and g is the gravitational constant. The term 6prgw describes the Stokes friction, where r is the particle radius, g is the molecular viscosity coefficient. The density of water, q water , comes directly from the physical ocean model while the density of the cod eggs, q egg , is computed from the equation of state based on salinity values, at which the eggs are assumed to be neutrally buoyant in the water column. The observations of cod egg density closest to our study area are off the coast of Helgeland in northern Norway, which indicate a salinity of 31.2 6 0.7 (Jung et al., 2012). However, using these values results in a vertical distribution where most of the eggs are confined to the surface layer, which is not consistent with observations (Lough et al., 1996). Since we know that the ocean reanalysis used as forcing has a saline bias in the NCC, a higher salinity is as- sumed for the eggs, 33.5 6 0.7. This yields a deeper and more realistic vertical distribution of the eggs. The radii of the eggs, r, are set to 0.675 6 0.025 mm, also based on Jung et al. (2012). The molecular viscosity, g, is computed as a function of temperature, salinity, and depth (Riley and Skirrow, 1975). The IBM describes the egg and yolk-sac larval stages of Atlantic cod using the temperature dependence for development rates as described in Daewel, Peck, et al. (2011). The equations for endogenous feeding (egg and yolk-sac larvae) include temperature-dependent utilization of yolk reserves and concomitant somatic growth of larvae. The effects of temperature on the duration of the egg stage and yolk-sac larval stage of Atlantic cod were taken from Geffen et al. (2006) and Jordaan and Kling (2003).

The model forcing Ocean reanalysis and derived fields
As the oceanic input to the oil and larvae drift simulation, the output from the TOPAZ reanalysis (Sakov et al., 2012) is used. This model system is based on the Hybrid Coordinate Ocean Model (HYCOM: Bleck, 2002). The model assimilates remotely sensed sea surface temperature, temperature, and salinity profiles from ARGO floats and sea ice concentration. The atmospheric forcing is the ERA-Interim (Dee et al., 2011) and data are assimilated using the Ensemble Kalman Filter (EnKF: Evensen, 2009) with 100 ensemble members (Sakov et al., 2012). The reanalysis fields provide daily information on currents, temperature, and salinity (Table 1), which are the ocean-fields required to force the oil drift model. The larvae drift model also requires information about vertical velocity and vertical diffusivity. Vertical diffusivity was not archived during the reanalysis production and vertical velocity is not computed in HYCOM, therefore both of these were diagnosed from the available model fields. The vertical velocity is estimated by integrating the continuity equations downwards in isopycnal coordinates as described in the HYCOM documentation (https://hycom.org/attachments/067_vertical_vel. pdf). The vertical diffusivity is diagnosed using a Richardson number based model (Pacanowski and Philander, 1981).
where k is the vertical eddy diffusivity, v the eddy viscosity, and Ri is the Richardson's number. The parameters chosen for the estimation of the diffusivity The temporal resolution of the ocean reanalysis output is daily and the spatial resolution is $12 km on a curvilinear grid. The model fields are projected from the native curvilinear model grid onto a regular 0.25 Â 0.125 longitude-latitude grid and nine depth-levels (Table 1) before being used as input to both the oil drift model and the IBM.

Wave forcing
The oil drift model needs the Stokes drift as well as the significant wave height and mean period. These variables were obtained from a global hindcast with WaveWatch3 with parameterizations described in Ardhuin et al.(2010) (The product is available for download at: ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST.). The temporal resolution of the wave input is 3 h and the spatial resolution is 0.5 (Table 2).

Atmospheric forcing
The wind field at 10 m from the ERA-Interim (Dee et al., 2011) is used to compute the direct wind-driven component of the oil drift model in OD3D. This is the same atmospheric fields were used to force the ocean reanalysis. The spatial resolution of the atmospheric fields is $80 km and the temporal resolution 6 h ( Table 3).

Spill locations
The oil drift model simulated oil spills at 49 different locations ( Figure 1), which are selected based on information about shipping routes in the study area. Automatic Identification Systembased information on ship density can be viewed on the webpage operated by the Norwegian Coastal Administration (https:// kart.kystverket.no/). Since only image data-not numerical data-were available the oil-release locations were determined by analysing the image data.

Spawning location and timing
As input to the IMB the time and location of the spawning of Atlantic cod in Norwegian waters were obtained from maps of spawning locations (downloaded from the Norwegian Directorate of Fisheries http:/kart.fiskeridir.no/). The dataset contains the geographical boundaries, as well as the starting and ending months of the spawning period for each spawning ground along the Norwegian coast. This information was remapped to the same grid as the oceanic forcing and the resulting fields ( Figure 2) show the locations where eggs are released in the model. The spawning intensity at the spawning locations is approximated as a Gaussian curve centered between the first and last month of spawning (Figure 2, right). In reality, the intensity of the spawning would vary between spawning grounds and from year to year (Sundby and Nakken, 2008; Opdal and Jørgensen,

The statistical risk calculation
The position and state of each of the released oil particles are recorded daily to compute the dispersion of oil from that point. For each spill location and month, the probability of finding oil in a certain location is computed on a 0.1 Â 0.1 degree grid. Since "aged," "evaporated," and "permanently submerged" oil is not tracked by the simulation, only "oil at the surface," "dispersed oil," and "stranded oil" is included in the probability calculation. "Stranded oil" would not affect the cod eggs or larvae, but was included in the probability calculation in case we would apply the risk maps to potential contamination of other species. A range of probability was decided upon to correspond to different levels of risk of oil-exposure: "low" is between 1 and 10%, "moderate" is between 10 and 30%, and "high" is more than 30% probability of the oil spill ending up in that particular region. The overlap between oil and eggs/larvae is computed by gridding the eggs and larvae onto the same 0.1 Â 0.1 degree grid that the probability is computed on. The indicator for contamination used in this study is the number of eggs and larvae that are present in a grid cell that contains oil. Depth is also taken into account when computing the overlap i.e. the eggs or larvae must be in the same depth interval as the oil in the model to be counted as contaminated. The depth intervals used are defined by the depth levels of the ocean model fields (Table 1). The method can differentiate between the eggs and larvae that are inside an area of high, medium, or low risk for oil exposure, but in the "Results" section we base the results only on the numbers of eggs and larvae located in regions with any amount of oil. To provide a relative measure of the risk, the total number of eggs or larvae that is in a region with oil spilled from one location in a certain month is divided by the average number of eggs or larvae overlapping with oil for all the oil spill locations in the period from January to May. This gives a measure of the relative risk for a certain spill location and month compared to the average risk and we will refer to this number as the "contamination risk."

Environmental variables
The contamination risk is related to the environmental variables that are used as forcing for the oil drift and IBM. The atmospheric and oceanic variables considered are listed in Table 4. Monthly time series of these variables were computed for each spill location by calculating the monthly mean condition over a 4 longitude and 2 latitude box surrounding each spill location. The geographical extent of the oil spill, larvae, and eggs (in km 2 ) as well as the total number of particles were also computed for oil, larvae, and eggs to determine if there is any connection between these properties and the overlap between oil and early life stages of cod.

Results
The largest concentrations of early life stage of cod occur in March and April, $79% out of the total number of eggs and larvae being present in these 2 months. Of the remaining 21%, $14% or the eggs are present in February and $13% of the larvae in May. Based on these numbers we focus on the analysis of the results in March and April, first presenting the environmental variables in these 2 months, then the general distributions and dispersion of eggs, larvae, and oil. In the "Risk of contamination from different spill locations" section, we present an analysis of which spill locations posed the greatest risk to early life stages of cod and finally analyse which environmental variables determined the variability in the risk of contamination in the "Physical factors determining the risk of overlap between oil and fish larvae" section.  The environmental variables that we consider in relation to the contamination risk are not all independent. It is therefore useful to understand the relationship between the different variables before considering the correlation with oil contamination. The Pearson's correlation coefficient is computed between the environmental variables used as forcing. The environmental variables that we consider are listed in Table 4 and the correlations between the variables shown in Figure 5. The highest correlation coefficients are between wind and Stokes drift and between wave height and wave period. Naturally, there is also a high correlation of the wind speed and direction with the ocean currents. Temperature and salinity are, surprisingly only correlated with the wave parameters, which are not used as an input to the ocean model. However the connection may be through a correlation between waves and westerly currents, which would bring warm, salty water into the coastal areas. There is a negative correlation between the meridional winds and the vertical velocities. This can be explained by coastal upwelling/downwelling; persistent winds from the south would cause downwelling, while winds from the north would result in upwelling along the west coast.

Distribution of eggs, larvae, and oil
To understand the results on contamination risk it is useful to understand the general temporal and spatial distribution of eggs and larvae. Though there is some interannual variation in the spatial distribution, the interannual variation in the total number of larvae is generally small. Around the time of maximum number of eggs/larvae one standard deviation is $10% of the mean number of eggs at that time and $20% of the mean for larvae. The model assumes that the spawning grounds, spawning effort and timing are the same each year which results in the seasonality of eggs and larvae abundance also being quite similar from year to year ( Figure 6). Moreover, the water temperatures from January to June, which affect the stage durations, are comparable between years. A typical temperature in this area is 7 C, at which the duration of the egg developmental stage is $13 days and the yolk-sac larval stage is $8 days. In an additional experiment where the currents are the same each year, but temperature and salinity vary between the years the standard deviation is only 1-2% of the total number of eggs or larvae. This is a clear indication that nearly all the inter-annual variability in the time series originates from the variability in the ocean currents and not from variability in temperature or salinity. The general spatial distributions of both eggs and larvae are also very similar in both March and April (Figure 7). The distribution is primarily decided by the currents and the total number of eggs and larvae, which are both similar in the 2 months. Although some currents are stronger in March, these currents are located further offshore than the regions with the highest concentrations of eggs and larvae.

Oil contamination of fish eggs and larvae
The modelled oil is released from discrete offshore locations. The areas of risk from each oil spill are elongated shapes around the spill site with more risk of encountering oil in the direction of the mean ocean current (Figure 8). In general, the difference between March and April is that for the spills along the west coast the spatial spread is larger and the high-risk areas are generally closer to the west coast in March (Figure 8). The stronger winds and Stokes drift in March are responsible for this. Oil is located mainly at the sea surface and hence influenced by the more variable environmental forcing such as winds and waves. Eggs, due to their positive buoyancy, are also generally located close to the surface and although they are only affected indirectly by the wind (through the ocean currents) the spatial variability of larvae/eggs is comparable to oil. The main difference is that the eggs and larvae are released in close proximity to the coast and hence can only be transported away from or parallel to the coast.

Risk of contamination from different spill locations
The location of the oil spill is an important factor for the contamination risk for both eggs and larvae (Figure 9). In general, the spills in the Skagerrak result in very low contamination risk, while the spills off the southwest coast of Norway have a very high probability of affecting the cod spawned in that area. Further north it is, not surprisingly, the spill locations closest to the coast and thus closest to the spawning ground that pose the highest risk. Overall, it is March that has the highest risk of overlap, although the total number of larvae is actually largest in April. For eggs the highest risk is in March, coinciding with the largest number of eggs. On average the contamination risk for eggs and larvae is $4 and 2.5 times higher in March than in April, respectively. When the total number of eggs or larva found in areas with oil is compared to the total number of eggs and larvae in the whole area it is clear that the eggs are about twice as likely to encounter oil, probably because they are on average located closer to the surface than the larvae.
In certain years and locations, the largest overlap occurs in April or February and for some individual spill location the amount of overlap varies dramatically from year to year. An example of the interannual variability in the risk of overlap between eggs and larvae and oil is shown in Figure 10. From this figure it is clear that an oil spill in a certain location can pose a great or small risk to impact the fish and larvae depending on the conditions in that particular year. In particular when both  the larvae and the oil are concentrated towards the coast, the resulting contamination risk is large. The interannual variability in the contamination risk is large for all the spill locations along the west coast of Norway ( Figure 11) and the contamination risk for larvae varies from almost 0 to 10. For some of the west coast oil spill locations close to the coast and in the northern part of the study area contamination risk is always greater than 1. At two locations in the southwest part of the domain, the lowest contamination risk is around 3. The interannual variability is larger in March than in April and for eggs, which have a lower abundance in April than in March, this difference is particularly pronounced. For the release locations in the Skagerrak the contamination risk is almost always zero, except for a very few years when oil spills from the area between Denmark and Norway have the same contamination risk as the highest values on the west coast.

Physical factors determining the risk of overlap between oil and fish larvae
To determine which physical drivers affect the overlap, the oil spill locations along the west coast are clustered into six groups based on the proximity of their location (Figure 12, left column). The oil spill locations in the Skagerrak are disregarded in this analysis since there is overall small contamination risk in this region ( Figure 9). The environmental time series are computed in a box surrounding each oil spill site as described in the "Environmental variables" section. Then the Pearson's correlation coefficients between the time series of the various environmental forcing factors and the monthly time series contamination risk are computed (Figures 12 and 13).
Wind affects the oil drift directly in the oil drift model. The results show that in groups 1-3 wind explains a similar amount of the contamination risk as Stokes drift. However, wind is highly correlated with both currents and the Stokes drift and to better understand the effect of wind an additional oil drift simulation where the wind is set to zero is performed. The contamination risk for this simulation is computed and the resulting correlations between environmental variables and contamination risk are estimated (not shown). In March, the run without wind has almost the same correlation coefficient as the ones shown in Figure 12, meaning that in March, it is the currents and Stokes drift that determine the variability of the contamination risk. In April, however, except for currents in group 3 there is no longer any environmental variable that explains more than 25% of the variability. Consequently, this shows that the wind forcing is central in explaining the variability in the contamination risk in April.
In March there is a correlation between the zonal current and groups 1, 3, and 6. This means that when there is an eastward current an oil spill in this area would tend to be transported towards the locations of the eggs and larvae. In April, the meridional currents appear to be more important and a slow, southward current would increase the overlap, while a fast, northward current would reduce the contamination risk. However much of this correlation is a result of the currents being correlated with westerly wind. This effect is present for most of the considered spill locations, but strongest for groups 1-3.
Salinity and temperature generally have small correlation with contamination risk. This is expected as most of the variability in the location and extent of both oil and eggs/larvae is determined by transport by the current and additionally Stokes drift and wind in the case of oil.
In March the westward Stokes drift has a correlation coefficient >0.6 with the contamination risk for both eggs and larvae in group 1 and >0.5 for eggs in group 2. Stokes drift is not included in the ocean forcing for the eggs and larvae, thus it is the transport of oil by the Stokes drift that can explain a large part of the variation in overlap. In April spills from groups 1 to 3 with simultaneous westward Stokes drift increase the contamination risk, but as mentioned before this is probably a result of the strong correlation with westerly winds. In the same month, a northward Stokes drift, just as the northward current, decrease the contamination risk, though the effect is weaker than for the currents.
The wave height and period are also correlated with the contamination risk. It is likely that there is no direct effect on the transport, but that the correlation is the result of the strong correlation between wave height and both wind and Stokes drift ( Figure 5). Finally, the vertical diffusion correlates positively in April with overlap between eggs and oil spills from groups 2 to 4. The vertical diffusion would tend to mix eggs deeper in the water column and this would affect their trajectory when there are vertical gradients in the currents. Typically, vertical gradients in the (baroclinic) currents will develop when the water column is stratified, which is more the case in April than in March.

Discussion
A model system is set up to investigate the risk of ship-born oil spill contamination in different locations and times for the early life stages of cod. The system consists of an oil drift model and IBM for cod early life stages that simulates larval drift as well as temperature-dependent development of cod eggs and larvae.   Figure 9. Contamination risk for eggs (left) and oil and larvae (right). The lines from the map on the right to the figure indicate the locations of the different spill locations.
exposed to oil from spills in a certain location and month relative to the overall risk from all oil spill locations and months in the spawning period.
The model results show that April and March are the months when the contamination risk is large, the primary factor being the high concentrations of eggs and larvae during these months. Oil spills close to the coast along the west coast pose the largest risk to the spawning grounds because of the proximity to the eggs and larvae. At the same time there is large interannual variability, depending on the conditions in individual years. A spill at one location can either not impact the eggs and larvae at all or contaminate a large area with high concentrations of eggs and larvae. Both the average contamination risk and the variability of that risk are much larger in March than in April, despite the fact that the total number of larvae is larger in April. March has on average stronger winds, waves, and currents and overall these results in more oil coming closer to the coast where the largest concentrations of eggs and larvae are found (Figures 8 and 9). The time series of monthly averages of the environmental forcing are related to the times series of contamination risk using simple linear correlation between the variables. For some oil spill locations there are correlations with some of the environmental variables that explain more than 25% of the variability. At those locations ocean currents and Stokes drift are the primary factors for the variability of the contamination risk in March, while in April the winds can explain the variability. At other locations there is no strong relationship between any of the monthly environmental variables and the contamination risk. One explanation for the lack of correlation with the environmental variables could be that the monthly averages do not retain the signal of important processes, for example the passing of a storm. The analysis of such highfrequency variability is left for future studies.
The principle interactions between oil spills and cod early life stages in different areas of the ocean has been shown in comparable approaches, combining oil drift, and IBMs (Spaulding et al., 1983;Vikebø et al., 2014). Vikebø et al. (2014), for example, included vertical migration behaviour of the cod larvae. Furthermore, their oil drift model provided concentration of polycyclic aromatic hydrocarbons and by following the larvae over time it was possible to record the integrated exposure to these chemicals. Those studies, however, only covered restricted time periods, and hence allowed neither for investigating the role of interannual variations in the environment for the impact of oil spills, nor for formulating a general risk index for oil spill contamination in cod early life stages. By estimating the overlap between oil spills and cod eggs and larvae for a 20-year time period, we can identify the role of interannually varying weather and ocean currents for the risk of an oil spill impacting cod early life stages. Although the oil spills simulated here are hypothetical, the approach creates a map of potential situations possible in an actual oil spill event. From our simulations, we can clearly identify both time (March and April) and location (close to the coast/ southern Norwegian Coast) where the risk for oil contamination is highest. We also found strong interannual variability in risk index that appeared mainly related to the wind and current situation, rather than to temperature-induced changes in growth rates of cod early life stages. The obtained knowledge on the correlation between environmental variability and overlap can then together with information on the actual situation in such an event be used to quickly estimate the risk for cod early life stages. However, to make this long simulation possible, certain model simplifications are necessary.
The spatial resolution of the environmental forcing is relatively coarse. This is, however, necessary to reduce computational demands for the long-term simulation especially for reanalysis runs that have the added computational burden of the data assimilation system. The atmospheric model runs with a resolution of $80 km in space, which resolves the wind associated with Oil contamination of fish eggs and larvae major weather patterns such as low-pressure systems. The temporal resolution of 6 h allows good representation of changes in winds associated with passing low-and high-pressure systems. The atmospheric forcing does not resolve finer-scale features associated with topography along the coast. The winds are used both directly to compute the wind-forced oil drift and to force the ocean model. No sensitivity study has been done for this model with respect to the resolution of the atmospheric forcing, but other studies on ocean models indicate that this affects the ocean currents (Jung et al., 2014). Changes in the wave condition were simulated with a resolution of 0.5 ($50 km), and with a temporal resolution of 3 h, which is sufficient to resolve the main changes in the large-scale wave conditions. The ocean model had a resolution of $12 km allowing the resolution of the main currents in the area, but some finer-scale features of the ocean's interaction with the coast are not resolved. In addition, this model permitted only larger eddies, while the smaller eddies that we know exists quite numerously along the Norwegian coast (Sundby, 1984) are not represented. The assimilation of sea-surface anomalies from altimetry into the model however corrects the large-scale currents (Sakov et al., 2012) and allows for a realistic representation of the ocean state. The ocean reanalysis should therefore captures the mean state of the system better than a free run with higher resolution. The daily resolution is sufficient to resolve the mean ocean currents, which changes more slowly than the atmospheric conditions. The exception to this is the tidal currents, which are not included in the present model simulation.
Considering the various resolutions of the forcing, the result that there is little-to-no risk of an oil spill in the Skagerrak impacting the early life stages of cod must be interpreted with caution because in the wave model the coastline is not well represented. In addition, our simulation has very few spawning grounds along the coast in the Skagerrak.
Tides are neglected in this study because it is challenging to assimilate sea-level anomaly in a model that includes tides, this is because in the observed sea-level anomaly that is assimilated the tides have been filtered out. Therefore, tides are not included in the ocean reanalysis used here. The tidal amplitude along the coast of South Norway is relatively small (<0.7 m with the highest amplitudes in the northern region) and so is not associated with strong currents, except in particular regions such as narrow fjord inlets, which are not resolved by our models. It is not possible to quantify how adding tides would change the outcome of this study and there is some uncertainty connected to the omission of the tides. Tidal currents vary over the day and in amplitude over the month, but the interannual variability is small, so excluding tides can affect the overall contamination risk, but is not likely to influence the interannual variability and its relation to the environmental forcing. Ocean surface waves induce two different processes in the upper ocean: Stokes drift (Stokes, 1847) and Langmuir circulation (Leibovich, 1983). The Langmuir circulation effects are not included in any of the models, they are observed to affect the spreading of oil (Simecek-Beatty and Lehr, 2017). The Stokes drift is included in the oil drift model applied here via the wave model forcing, but not in the forcing for the IBM for the cod eggs and larvae. However, it has been shown that both of these waveinduced processes affect the drift of fish eggs and larvae close to the surface (Röhrs et al., 2014). Recently these effects of ocean waves are also included in general circulation models (GCM) (Staneva et al., 2017) including the ocean model used here, HYCOM (Breivik et al., 2016), but are not yet available as part of ocean reanalysis. It will probably be some time before wave processes get included in ocean reanalysis products, but future studies can potentially include wave effects also for the egg and larval IBM.
Some uncertainty is also related to the spatial resolution applied in the analysis methods. When computing the probability maps for oil and the overlap in locations among oil, eggs, and larvae, all variables were mapped onto a 0.1 Â 0.1 grid. This was chosen as it is reasonably close to the ocean model resolution, which is the forcing with the highest spatial resolution. A test was done for computing risk maps using a finer grid (0.05 ) and gave similar results, so 0.1 was chosen to be adequate. The criteria used for a grid cell to be counted as "affected by oil" an egg or larva needs to be in the same 0.1 Â 0.1 grid cell and at the same depth layer as an oil particle. Obviously fish eggs and larvae as well as the oil particles are much smaller than the size of these grid cells, so being in the same cell does not necessarily mean that the eggs or larvae actually interact with the oil. However, in the Lagrangian-based modelling approach individual particles typically represent super-individuals, which encompass a multitude of the particles actually simulated (Grimm and Railsback, 2005). Since this is also true for the oil particles, the assumption that the simulated particles overlap in the considered   Correlations that explain more than 25% of the variability are marked with black circles, 'þ' marks positive correlation.
grid cell is still valid. In addition, with the uncertainty related to the coarse resolution of the ocean current and wind fields, and the model error in general, a stricter criterion for overlap is not necessary, and would require a larger number of simulated particles to obtain robust results. The overlap between eggs and larvae varied quite dramatically between years even when the interannual variability of spawning itself is neglected. Sundby and Nakken (2008) found a relationship between roe-index, an indication of spawning intensity and the temperature in the Kola section in the Barents Sea. However, only one of the spawning grounds they considered is within our study area and the study period was prior to 1980. Other studies, also from the period before our simulation period have attributed the spawning variability to size-selective fisheries (Opdal and Jørgensen, 2015). None of the studies offer any information about interannual variability in the timing of the spawning, thus we did not have sufficient information about the spawning variability for consideration in our analysis. Better knowledge of what causes interannual variability in the spawning would improve the estimates of risk from our study. However, keeping the spawning constant between years also makes it easier to isolate the role of the environmental forcing for the interannual variability.
Recently, Langangen et al. (2017) showed that the impact of oil spills on the cod recruitment is additionally related to spatial differences in cod natural mortality, an effect that also should be considered in future studies on oil spill impact on fish. Many marine animals exhibit seasonal behaviour, and the oil drift also takes different routes depending on the season. Studies such as this can therefore be used to determine in advance whether certain species or their reproduction are endangered by pollutants in the ocean. For species that are not drifting with the currents, but with known migration patterns, the oil-risk maps can be used directly without the need for an IBM. Our study provides an improved understanding of how oil spills from certain locations can be a risk to early life stages of cod in different seasons and under different environmental conditions based on available ocean and atmospheric forcing. Both the quality and resolution of reanalysis products are constantly being improved as more computational resources become available. Satellite and autonomous instruments also collect more observations that can be assimilated and used for evaluating and improving these models. Thus, the outlook for developing improved estimation of the drift of both organisms and pollutants, and getting a deeper understanding for the underlying mechanisms behind what brings them in contact with each other is very good.  Latitude Figure 13. Correlation coefficient between 20 years of contamination risk and the respective environmental variables in April. Left column: Geographical location of the release points that are used to relate the contamination risk to the environmental variables. Right column: correlation coefficient between contamination risk for eggs and larva and the environmental variables in the area of the of oil spill locations.
Correlations that explain more than 25% of the variability are marked with black circles, 'þ' marks positive correlation.