Quantifying the consequence of applying conservative assumptions in the assessment of oil spill effects on polar cod (Boreogadus saida) populations

In order to assess the potential impact from oil spills and decide the optimal response actions, prediction of population level effects of key resources is crucial. These assessments are usually based on acute toxicity data combined with precautionary assumptions because chronic data are often lacking. To better understand the consequences of applying precautionary approaches, two approaches for assessing population level effects on the Arctic keystone species polar cod (Boreogadus saida) were compared: a precautionary approach, where all exposed individuals die when exposed above a defined threshold concentration, and a refined (full-dose-response) approach. A matrix model was used to assess the population recovery duration of scenarios with various but constant exposure concentrations, durations and temperatures. The difference between the two approaches was largest for exposures with relatively low concentrations and short durations. Here, the recovery duration for the refined approach was less than eight times that found for the precautionary approach. Quantifying these differences helps to understand the consequences of precautionary assumptions applied to environmental risk assessment used in oil spill response decision making and it can feed into the discussion about the need for more chronic toxicity testing. An elasticity analysis of our model identified embryo and larval survival as crucial processes in the life cycle of polar cod and the impact assessment of oil spills on its population.


Introduction
The exploration, production and transport of oil and gas in the Arctic region involves challenges related to environmental, socio-economic and cultural issues as well as good governance. High quality knowledge about the potential impact of exploration, production and transport activities and adequate use of this knowledge is a prerequisite to make informed decisions in case of an oil spill. An important instrument within this context is the net environmental benefit analysis (NEBA). This is a procedure for assessing the environmental consequences of oil spill response options. It enables decision-makers to select appropriate responses minimizing potential impact (Wenning et al. 2018). NEBA relies on sound risk assessments using oil toxicity data that assessed potential oil impacts on relevant biological resources preferably at the population level. A worst-case (precautionary) assumption is often applied within a NEBA using fixed threshold values. These threshold values cover both short and long term effects of oil exposure (e.g. Coelho et al. 2015), although in practice acute data usually form the main basis for these threshold values. The current study aims to explore the consequences of using these threshold values as cut off points by assessing population level impacts of oil exposure for polar cod (Boreogadus saida) in the Arctic. Polar cod is the most widespread and abundant fish of the Arctic Ocean (Lowry and Frost 1981;Parker-Stetter et al. 2011) and it plays a key role in the Arctic marine ecosystem (e.g. Mueter et al. 2016) connecting zooplankton with marine vertebrate predators, such as marine mammals and birds (Welch et al. 1992). An oil spill effect on B. saida populations will therefore potentially cascade through the Arctic food web, as demonstrated for sea ice effects where modelling results suggested negative top-down effects by polar cod on amphipods (Stige et al. 2019).
Polar cod is well-adapted to large annual fluctuations in physical (ice cover) and biological (primary productivity) characteristics, which probably explains its overwhelming success in Arctic waters (Lowry and Frost 1981). Due to its ecological role and abundance, it is a good model species to study.
In the current study, the toxic effect of oil was investigated with two different approaches: a precautionary approach, in which all exposed individuals survive when exposed to concentrations below a threshold value and die when exposed above this threshold value, and a refined (fulldose-response) approach, in which an assumed relationship results in an effect per toxic unit (TU). An age-based matrix model was developed to simulate the population development of polar cod. A benefit of matrix models is that they are analytically tractable and their behaviour is well analysed (Caswell 2001). They are often used in ecotoxicology and risk assessment approaches (e.g. Caswell 1996; Klok and De Roos 1996;Smit et al. 2006;Klok 2008;Bergek et al. 2012;De Vries et al. 2018).
With a daily-based simulation of the matrix model we investigated the impact of selecting either the precautionary approach using a single threshold value or a refined approach using a full dose-response relationship on the population dynamics of polar cod. A similar exercise, published by De Vries et al. (2018), was performed for another well-adapted, keystone organism in the Arctic Ocean, namely for herbivorous copepods of the genus Calanus (Welch et al. 1992;Falk-Petersen et al. 2007).

Overview of steps
To investigate population level effects of oil spills on polar cod in the Arctic, an age-structured matrix population model was developed. In this model, the life-history parameters survival and reproduction were parameterized with values and relations taken from literature. As an important step in the modelling process we performed an elasticity analysis for the underlying parameters to assess the effect of a relative change in each of the life-history parameters on the relative change of the population growth factor. To incorporate toxic effects of oil into the population model, mortality and/or reproductive output endpoints are required. Because the relationship between exposure and the effect on the life history parameters of polar cod was not available in literature, we used hypothetical effect levels based on an acute to chronic ratio (ACR) (for fish) in this study to scale up to the population level (May et al. 2016).
The (scarcely) available literature on toxic effects of oil on polar cod vital rates was only used to assess their relation to our hypothetical effect levels, by calculating field relevant toxic units based on these data.

Life cycle and model approach
For this study we used a matrix population model as these models: (1) have a direct relation to empirical age-structured data, (2) the resulting population growth factor is clearly linked to life-history parameters such as reproduction and survival, and (3) a relatively low number of input data is required (Beissinger and Westphal 1998).
First, literature was surveyed to establish the life cycle of B. saida (reported in detail in Online Resource 1). The life strategy of polar cod consists of early maturity, relatively rapid growth in comparison to other gadoids at cold temperatures (Laurel et al. 2016), small body size, and a relative short life span that rarely exceeds 5 years with a maximum of seven (Craig et al. 1982;Cohen et al. 1990;Jensen et al. 1991). It can be found in demersal, pelagic, and sympagic (in association with sea-ice) habitats, depending on the life stage (Cohen et al. 1990). Aggregation of individuals (schools) is known (Welch et al. 1992;Crawford and Jorgenson 1996). Figure 1 represents the life cycle on which the model was based. It includes the assumption that the male to female ratio equals 1:1 and only the female part of the population was considered (by basing the model on female vital rates), thereby discarding the female-skewed sex ratio in later years. Also a region with optimal conditions was assumed enabling polar cod to reach maturity in 3 years and a maximum age of approximately 7 years. For simplicity, constant yearly spawning was assigned to take place in a single day for fish in their fourth to seventh year.
This results in a post-spawning population model with a generation time of 4 years (x 4 (t) is the class that represents this population part). The four classes (x 1 (t), x 2 (t), x 3 (t), x 4 (t))) T represented newly hatched polar cod (0+), 1 year old, 2 years old and older than 3 years. The Leslie matrix M for this system was given in Eq. 1. Time t was measured in years. Further details are presented in Online Resource 2 ( Fig. S2.2) such as the structure of the year and how this is implemented in the model, the derivation of the matrix elements from daily based survival probabilities, and what these elements look like in terms of daily survival and fecundity (m). The yearly survival probability in class x i (i = 1, 2, 3, 4) is y i and y i = s 365 i with s i the daily survival probability. Daily embryo and larval survival was denoted as s 0 ; b denotes the hatching duration of the embryos in days.

Parameterisation
We based model parameter values (Table 1) on information obtained from literature making several assumptions. The temperature tolerance of polar cod adults is much higher than larvae. Adults are found in waters with temperatures ranging from −2 to 13.5 °C . Egg and larval development of polar cod require temperatures below 3 °C for normal development of eggs and larvae (Aronovich et al. 1975;Sakurai et al. 1998;Dahlke et al. 2018) and is under pressure from global warming (Dahlke et al. 2018). Nahrgang et al. (2015) suggested that polar cods are total spawners (i.e. release one batch of eggs per spawning season) with a group-synchronous development. The duration of the egg laying period is therefore assumed to be 1 day (see section 'Reproduction' in Online Resource 1). The number of eggs produced per female increases with age and totals 75,000 eggs per female over their entire life span . Most females generally mature at age 3 (Craig et al. 1982;Cohen et al. 1990). For the model, the egg production was assumed to be constant over the years and reproduction took place from the fourth until the seventh year of life. This resulted in an egg production of 18,750 eggs per female per year. The daily embryo and larval survival was based on data from Sakurai et al. (1998) and Nahrgang et al. (2016), and was on average 1% (see also Table 2). Hatching duration (b) also depends on water temperature. Data from Sakurai et al. (1998) were used to describe the temperature dependent hatching duration (see Online Resource 2). The daily survival in the juvenile (including older larvae), sub-adult and adult stages of polar cod were calculated using the relationship between body weight and natural mortality in juvenile and adult fish (Lorenzen 1996) (see section 'Mortality' in Online Resource 1). After assigning the parameters a value based on the literature an elasticity analysis  Aronovich et al. (1975) and Melnikova and Chernova (2013) b Based on data from Sakurai et al. (1998) c Based on data from Nahrgang et al. (2016) d Based on the weight dependent mortality relationship (Lorenzen 1996, Online Resource 1, Eq. S1.1) d Based on Hop et al. (1995) and Hop and Gjøsaeter (2013)  The number of eggs laid in one day per female 9 × 10 3 to 6 × 10 6d,e 18,750 e No was performed to determine the relative sensitivity of the model to changes in parameter values (Online Resource 2).

Toxicity
Toxic units (TU) were used in this study to express oil exposure (Von der Ohe and De Zwart 2013). To derive exposure values x (in TUs, see details in Online Resource 3) from a (measured) oil compound concentration c was scaled with its inherent effect concentration (e.g., acute 50% lethal concentration, LC50 4 , for adult fish (= class 4)): x = c∕LC50 4 . Both the No Observed Effect Concentrations (NOEC) and 50% lethal concentrations (LC50 i ) are also scaled with the LC50 4 (for later use in Eqs. 2 and 3): n = NOEC/LC50 4 = 1∕ACR , L i = LC50 i ∕LC50 4 for (i = 0, 1, 2, 3, 4, where 0 represents embryo and young larvae up to the time that all embryos are hatched). Thus L i is the LC50 in stage i relative to the LC50 in the adult stage. A theoretical relation based on the acute-to-chronic ratio (ACR) for fish (May et al. 2016) was used to parametrize the concentration-time-response-relationship on the population growth factor. These TUs were used to express the exposure to both single oil components and to mixtures of oil components.
Although oil toxicity could affect several of the life history parameters we focused on toxic effects on survival. Because limited toxicity data were available, a theoretical approach was used to describe the relation between exposure and effect. We used the same approach as described by De Vries et al. (2018) assuming that the hazard rate (h(t), the probability per unit of time to die at time t), is affected by an exposure x (in TU) above the No Observable Effect Concentration (NOEC) of n TU: with the h 0 (t) baseline or natural hazard rate (including predation, see Lorenzen 1996). The relationship is assumed to be multiplicative, and the magnitude of the effect is expressed as exp(β i ) per toxic unit above the NOEC. Following De Vries et al. (2018), the effect magnitude is expressed as (see details in Online Resource 4): A difference with the approach of De Vries et al. (2018) is that the LC50 (L i , expressed as toxic units) was assumed not to be equal in each life stage. Ratios of effect concentrations between fish life stages as reported by Hutchinson et al. (1998) were used for this purpose (Table 3). An ACR value of 12 was used for fish (as derived by May et al. 2016) and an acute exposure time t a of 4 days (as per OECD standards for fish tests (OECD 1992)).
Crude oil is a complex mixture of both hydrocarbons, such as alkanes, cycloalkanes and aromatic hydrocarbons, and non-hydrocarbon compounds, with varying effects on exposed biota. Toxicological risks of oil mixtures are primarily determined by its water soluble components (e.g. French McCay 2002; Olsen et al. 2013). Information on the impact of oil-components on life-history characteristics in different developmental stages are relatively scarce (Olsen et al. 2013). Some effect values on vital rates are available for polar cod exposed to specific crude oils and were only used in the discussion to place the modelling results into context of field relevant concentrations (reviewed by De Hoop et al. 2011;Gardiner et al. 2013;Nahrgang et al. 2016) (see Online Resource 3, Table S3.1).

Oil spill scenarios
For the model we consider a population of polar cod to be all individuals that are present in a specific region during the simulation/observation time. This population is considered to be a closed system. The model simulated the polar cod population on a day to day basis (Online Resource 2) exposed to a range of TUs for a range of exposure durations (2, 4, 8 and Total length (L in mm) of 1 to 4 year old polar cod in the Arctic region (mainly North-eastern part of Svalbard) , calculated wet weight (W, using Online Resource 1 Eq. S1.2), calculated mortality (using Online Resource 1 Eq. S1.1) and calculated daily survival of (sub)adults (based on the average mortality of 1 and 2 year olds (sub-adults) and the average mortality of 3 to 7 year olds (adults))  16 days) and a range of exposed fractions of the polar cod population. The exposure concentration varied between 0.9/ ACR (i.e., 90% of the NOEC) and 1.1 TU (i.e., 10% above the LC50) in 16 equidistant exponential steps and the exposed fraction of the population varied between 1% and 99% in nine logarithmic steps. The affected fraction indicated the fraction of the population that is being exposed to oil. Only this fraction was affected, following one of two different approaches evaluated. The approaches were: (1) precautionary: all individuals exposed above the NOEC die instantaneously; (2) refined (fulldose-response): individuals die as the result of an increased hazard rate, which depends on the exposure concentration and duration (as described above). Because the model is not spatially explicit, no exchange of individuals between the exposed and the unexposed fraction of the population was modelled. Recovery time was evaluated for the entire adult population consisting of the unexposed individuals and those that survived exposure.
The number of adult individuals peaked each year at the census moment (after hatching). The recovery time was defined as the time it takes to reach at least the number of adults of an unexposed population measured at the census moment after an exposure. Because this definition is linked to the census moment, which takes place just after hatching, the recovery time was always expressed in number of full years starting from year one of the simulation. Obviously, in reality, recovery can take less than a year.
For both the precautionary and the refined approach, results were obtained through simulation runs with the matrix model as given in Online Resource 2 (and Eq. 1). For this purpose, a period of 20 years was simulated, with oil exposure taking place during the breeding season at the end of the second simulated year. We report about the simulation where all life stages (including the embryo-stage) were affected by the exposure. We also considered a case where the first class (newly hatched polar cod 0+) was the only year class exposed to the oil spill and a case where all but the reproducing adults were exposed to the oil spill (the results of these last two cases can be found in Online Resource 5).
For each combination of exposed fraction, exposure concentration, and exposure duration, the resulting recovery times were compared for both approaches (precautionary and refined) and expressed as the ratio R of recovery time for the precautionary approach 1 (T r1 ) to the recovery time for the refined approach 2 (T r2 ):

Elasticities
The population growth factor based on the parameter values as used in the model ranged from 1.92 to 2.12 with temperatures varying between −2 and 3 °C (Fig. 2a). Obviously, the growth factor depends on the temperature because hatching time depends on temperature (see Fig. S2.1). The elasticities (see Fig. 2b, c, d and Online Resource 2) showed that the model, in terms of the population growth factor, was most sensitive to the daily survival of the embryos and larvae, and less to the yearly survival in the fourth year class, whereas the variation in length of the hatching period b had only a minor impact on the population growth. The production of embryos and the survival in the first years of life had the same elasticities (Online Resource 2). This suggests that the embryo and larval survival are crucial factors to consider in oil toxicity experiments. Figure 3 shows the ratio of recovery times of the precautionary and the refined approach. A higher ratio R represents a larger difference between the precautionary approach and the refined approach with a more optimistic estimate for the refined approach as R>1. Figure 3 furthermore shows that temperature has only a minor effect on this relative difference. At most there is 1 unit difference in R between the simulated temperatures (ranging from −2 up to 3 °C).

Precautionary versus refined approach
We assumed that no effects occurred below the chronic NOEC for both the precautionary and the refined approach, explaining why the ratio of recovery times is by definition 1 for exposures below the NOEC in Fig. 3. For low concentrations (just above the NOEC of 1/12 TU), the difference between the precautionary approach and the refined approach were the largest (Fig. 3). For higher concentrations, the difference between the model simulations of the refined and precautionary approach decreased (Fig. 3). This was also the case for longer exposure durations (Fig. 3).
When less than half of the population is exposed (<51%), no difference in recovery time between the two approaches was detected with the model used here. When more than half of the population (≥51%) is exposed, the precautionary approach resulted in a longer predicted recovery time. The ratio of recovery times between the precautionary and the refined approach was eight at maximum (i.e. for the 99% exposed). This means that for these cases the recovery time in the precautionary approach was eight times longer compared to the refined approach.
We have also determined what happens if not all year classes are exposed to the oil spill (see Online Resource 5). If only the youngest age class is exposed to the oil spill, then the precautionary and refined approaches give the same recovery durations (so R equals one for all doses and exposure durations). If all age classes except the reproducing class are exposed to the oil spill then we see a similar pattern in the ratio of the recovery durations (R) as in our scenario where all stages are exposed to the oil spill.

Discussion
Our model for polar cod represented a population with an approximate doubling of the population every year at temperatures between −2 and 3 °C. This implies that the population is buffered by loss due to predation from higher trophic levels in the Arctic food web or other factors. The survival of embryos and larvae during hatching appeared to be of crucial importance for the population growth factor, as revealed by the elasticity analysis. With respect to the comparison of the precautionary and refined scenarios this study showed that the refined approach provided shorter recovery durations compared to the precautionary approach after exposure of all year classes to an oil spill. If only the youngest class is exposed then both approaches give a similar result.

Limitations of our study
The ability of polar cod to adapt to local conditions implies that its biological characteristics (biomass, feeding, growth rates, spawning and hatching) vary geographically, seasonally and from year to year (Lowry and Frost 1981;Jensen et al. 1991;Hop et al. 1995;Bouchard and Fortier 2008). The current understanding of polar cod is surprisingly fragmented and inconclusive, leaving major gaps in knowledge that prevent a holistic understanding of the interaction between the species and its environment (Mueter et al. 2016).
Our literature study (see Online Resource 1) revealed that (1) the spawning strategy of polar cod has yet to be fully understood (Nahrgang et al. 2015), (2) the mortality of embryos and larvae reported in literature (Sakurai et al. 1998;Fortier et al. 2006;Nahrgang et al. 2016) vary considerably, (3) the natural mortality of (sub)-adult polar cod is not measured, and (4) there is insufficient experimental information on the relation between chronic oil exposure levels and effects on survival in polar cod.
Based on earlier findings (Nahrgang et al. 2015) we modelled the yearly spawning (point (1) above) as occurring on one day, namely releasing one batch of eggs per spawning season with a group-synchronous development. To assess the effect of this assumption we have simulated also a population with spawning occurring over 10 days. This did not Fig. 2 The temperaturedependent model results. The model gives as a consequence of temperature dependent hatching period, at different ambient temperatures of the sea water differences in a the dominant eigenvalue, b the elasticity of hatching period b, c the elasticity of the daily egg and larval survival ( s 0 ), and d the elasticity of the yearly adult survival (note that from the formulas S2.9, S2.11-2.13 in Online Resource 2 it is clear that e(y 1 ) = e(y 2 ) = e(y 3 ) = e(m)) affect our results in a qualitative way, although the maximum of the ratio R was a bit higher (11 instead of 8; results not shown).
To overcome points (2) and (3) we used the weight dependent mortality relationship for fish (Lorenzen 1996;Peterson and Wroblewski 1984) to calculate yearly survival. In Online resource 1 we compared this value to other estimates for this species and similar species.
With respect to point (4), the lack of experimental data on mortality due to oil (specifically LC50 values), we extrapolated LC50 values from experimental NOEC values in polar cod using the acute-to-chronic ratio (ACR) of 12.0 (May et al. 2016) (see Online Resource 3 Table S3.1). Effects of oil components on polar cod are available (see Online Resource 1) but these have limited value for use in our population model, as most studies used biomarkers (genes, enzymes, metabolites, DNA damage) as endpoints of which the quantitative relation to the vital rate survival is unknown (Nahrgang et al. 2009;2010a, b, c;Fig. 3 The relation between exposure concentration and the ratio between recovery times of the two considered scenarios. The ratio R (Eq. 4) equals the recovery time of the precautionary approach (=worst-case scenario) divided by the recovery time of the refined approach. It is shown as dependent on the exposure concentration for three water temperatures: −2 °C (a), 0.5 °C (b) and 3 °C (c). In nine panels the percentage of the population that is exposed is increased logarithmically from 1 to 99%, see header of each panel. The ratio R is also dependent on the duration of the exposure (denoted with lines with different markers). The left dashed vertical line resembles the NOEC (n = 1/ACR in TU) and the right dashed vertical line the LC50 for adult fish (L 4 = 1 in TU). The minimum recovery duration is one year. The ratio R is equal to 1 when the precautionary approach and the refined approach result in the same recovery duration, when R is greater than 1, recovery times calculated with the precautionary approach are longer than those calculated with the refined approach. The curves are not smooth, because the recovery durations (used to calculate R) are expressed as full years. At lower concentrations, lines corresponding to different exposure durations overlap, and only the line corresponding with the longest exposure duration is visible Christiansen et al. 2010;Geraudie et al. 2014;Bender 2015;Bender et al. 2016).
Obviously the extrapolated LC50 data we use is less informative than experimental species derived LC50 data. Such more realistic LC50 data, however, is not expected to affect our model results in a qualitative way, since the highest sensitivity of the population growth factor on daily survival of embryos and larvae is independent of toxicity data and both the precautionary and the refined approach are based on the same extrapolated LC50 data.
Given that the population growth factor was most sensitive to the daily survival of the embryos and larvae, we recommend to invest in relevant experiments in early live stages of polar cod. These experiments should focus on chronic exposure (duration at least 4 days, to allow steady state in toxicokinetics), and report the number of survivors over time as well as the LC50.
The lack of polar cod specific chronic LC50 data, combined with information on the number of survivors over time, deters a full toxicokinetics toxicodynamic (TKTD) approach such as the Dynamic Energy Budget approach exemplified in Atlantic cod Gadus morhua (Klok et al. 2014). As noted by Klok et al. (2014) the mechanistic DEB approach can help to reduce experimental effort, since DEB allows extrapolation of toxicity parameters out of the range of tested exposure durations, and allows extrapolation from one tested petroleum substance to others with the same mode of action (e.g. non polar narcotic). Moreover, DEB models have the added value that sublethal toxic effects such as reduction in growth of individuals can be extrapolated to the population level (Kooijman 2010;Klok and De Roos 1996). This seems particularly relevant for polar cod since available laboratory data already show reduction in length (which implies reduced growth over time) of larvae at NOEC concentrations  ; Table S3.1).
Given limitation in available toxicant-induced mortality data, our current approach, using the ACR for fish, to parameterise the concentration-time-response-relationship and express exposure in TU (von der Ohe and De Zwart 2013), proved useful. A similar approach was also used by De Vries et al. (2018) to investigate the consequences of assumptions in oil spill assessment on population impact of Calanus species. For fish, the LC50 value, based on parent naphthalene concentration, was 30 µg L −1 for juvenile polar cod exposed to the water accommodated fraction of oil (Gardiner et al. 2013). Comparing this LC50 value with NOEC levels on larval survival (37 days exposure to a crude oil water-soluble fraction with a sum of 26 polycyclic aromatic hydrocarbons 2.18 μg L −1 )), would result in a ratio of 13.8. Because this ACR is close to the applied generic ACR of 12, our theoretical approach is supported. In addition, ratios of effect concentrations between fish life stages (Hutchinson et al. 1998) were considered acceptable for refinement of our hypothetical effect levels, although these ratios do not specifically apply to oil components.
It should also be noted that we did not account for the longer lifespan of females in our female-based model which can be inferred from the sex ratio of polar cod that approaches 100% females in cohorts with higher ages. This was not possible because the generic relation between weight and natural mortality (Online Resource 1, Lorenzen 1996) does not consider intraspecific and interspecific variation, and thus no differences between sexes.
An example of differences between sexes is that female fish have generally a lower ethoxyresorufin-O-deethylase (EROD) activity than male fish during spawning (Larsen et al. 1992;Arukwe and Goksøyr 1997). Because EROD activity points to potentially metabolizing toxic compounds (Sarkar et al. 2006), female fish might be more sensitive to contaminants than males. Luckily, the elasticity of the survival parameters of juveniles and adults is low in our model, but we might have underestimated the recovery times in the refined approach due to this. Although a linear model as we now have developed is not accounting for density dependence, the most obvious effect of compensatory growth in young fish is that they might grow faster and become mature earlier. If it was included in the model it would have occurred under both considered scenarios and the ratio R of the return times would not have been affected.

Reality check
The TU range applied in this study increases up to 1.1, which corresponded to an exposure concentration slightly above the LC50. Exposure above LC50 values can be realistic in field situations, especially directly after and/ or near the source of a spill (Table 4). TUs in field situations based on concentrations from actual, experimental and modelled spills (see Online Resource 3 for derivation of TU) ranged between <0.00003 and 19.11, with most values below 0.1 TU. Whereas the concentrations in the field usually form a gradient, here we only simulated a constant exposure concentration for each scenario. According to our modelling results this meant that for field situations, the predicted impact based on either approach (precautionary versus refined) would in most cases lead to comparable results in terms of recovery times. However, the precautionary approach results in impact estimates that are up to eight times as severe as those determined with the refined approach. This is especially the case when a larger part (more than half) of the population is exposed and when temperatures are relatively low. Hu and Wroblewski (2009) developed an age-structured Leslie matrix model of the larger Atlantic cod (Gadus morhua) and used three different values (0.2, 0.3, and 0.4 year −1 ) for natural mortality resulting in age-specific probabilities of survival for cod (ages ≥1 year) of 0.819, 0.741, and 0.670, respectively. Indeed, these yearly survival probabilities are much higher than those assumed for polar cod (and used within this study). Hu and Wroblewski (2009) acknowledged that, based on weight dependent mortality, theoretically, the natural mortality for northern Atlantic cod older than 1 year should be approximately 0.2 year −1 . This value is similar to an estimated mortality based on the longevity relationship (following Hoenig 1983, natural mortality of G. morhua is 1.83 year −1 , see Online Resource 1). The mean annual natural mortality rate estimated by Aanes et al. (2007) using an age-structured model based on catch-at-age data and abundance indices, ranged from <0.2 to >0.7 with a mean value of 0.36. This modelled value is higher than estimated using the longevity relationship. Aanes et al. (2007) report that the estimates are rather imprecise. The component accounting for the majority of the fluctuations in natural mortality was caused by year-to-year variability (Aanes et al. 2007). Bogstad et al. (1994) studied cannibalism using stomach content data and virtual population analysis for Atlantic cod (G. morhua). They found an average yearly mortality due to cannibalism and competition of 0.32-0.64 year −1 averaged over 2.33 years. Hu and Wroblewski (2009) also included higher rates to investigate the effect of high natural mortality (possibly due to predation) on cod population dynamics and found that the population growth factor for G. morhua is more affected by survival in all ages than by the fertility values. The analysis of the current model study showed that the same applies for polar cod, with the highest elasticity for the survival of embryos and larvae.

Comparison with model approaches for Atlantic cod
Of course commercial fishing also affects the population growth factor and therewith the population's recovery time after an oil spill (Reed et al. 1984). Therefore, effects of fishing and potential oil spills should preferably be evaluated jointly using the same basic principles (Carroll et al. 2018). Carroll et al (2018) simulated oil spills in the Atlantic cod (G. morhua) spawning grounds to assess potential impacts of oil spills on the cod fishery. Reductions in survival for pelagic stages of cod were up to a maximum of 43%, resulting in a decrease in adult cod biomass up to a maximum of 12%. The study suggests that the diverse age distribution helps protect the adult cod population (Carroll et al. 2018).
For polar cod, commercial fishing has been reported to peak in 1971 at 348 thousand tonnes, but after 1984 world catches declined steadily (FAO 2017). The total stock size of polar cod in that period is unknown, but it has been measured acoustically since 1986 and has fluctuated between 0.1 and 1.9 million tonnes from 1986(ICES 2014. In that period the greatest increase in abundance has been observed between 2003 (0.3 million tonnes) and 2004 (1.1 million tonnes), and it further increased up to 1.9 million tonnes in (McBride et al. 2016. The FAO reports that polar cod populations are little affected by fishing because r-selected species can support high levels of fishing mortality and have a short recovery time (FAO 2017).

Recovery times
The recovery times in our model were based on population growth rates under ideal environmental and ecological conditions. In reality, actual recovery times could be longer, e.g. due to higher actual mortality (such as high predation, fishing or extreme environmental conditions). Our results included the effect of temperature on the population recovery time by addressing the temperature dependence in embryo and larval development. A relatively high temperature corresponds to a relatively short hatching duration (i.e. larval development is temperature dependent (Sakurai et al. 1998)), which, in our model, resulted in a high hatching success due to the cumulative daily survival probability and consequent increase in population growth factor. These direct effects are relevant when considering short term, local variance of environmental conditions, but when discussing the effects of temperature in relation to climate change there are other aspects to consider. A warmer climate will lead to changes in phenotypic traits including earlier maturation, smaller size, increased investment in reproduction at early age and in sum a reduced fecundity , and to changes in the ecosystem (species presence, interactions, habitat and timing and quality of seasonal cycles). It could be argued that a short-term temperature increase up to 3 ˚C would enhance polar cod population growth (as indicated in this study based on the temperature-dependent hatching time) whereas prolonged exposure to elevated temperatures (i.e. global warming) will ultimately reduce population growth due to ecological cascading effects (as suggested by Nahrgang et al. 2014).

Consequences for oil spill response
Differences between the precautionary and refined approach for polar cod were not as large as found in an earlier study with Calanus (De Vries et al. 2018). Because not all exposed individuals die in the refined approach the recovery times are smaller than in the precautionary approach. That is the reason why the difference between the two approaches is largest for exposures of short duration and low toxic concentrations. The fact that the refined and precautionary approach give the same recovery times for cases where less than half of the population was exposed show that the precautionary approach is a good starting point for handling oil spills that do not potentially affect a large part of a population. Using either the precautionary or refined approach has consequences for assessing oil spill response options (i.e., NEBA). When considering multiple mitigating measures for oil spill response, multiple environmental compartments are assessed (water column, water surface and shore). Comparing results for these different compartments can only be sensible when similar approaches have been used. Only using the same level of conservatism for each compartment leads to results that can be compared. From a precautionary and pragmatic point of view, it is best to apply the precautionary approach, as it offers sufficient conservatism and requires less data. However, when weighing spill response options, a more realistic impact assessment can be advised.
For modelling a realistic impact assessment more information is required on temperature-dependent reproduction and survival of differently aged specimens of polar cod. With respect to toxicity tests in the laboratory additional knowledge is needed on how different life stages and sexes are affected by crude oil mixtures and different chemical substances. It is also important to be able to link external and internal concentrations in experimental settings. That way oil spill scenarios can focus on external concentrations.