Impacts of Oil Spills on Arctic Marine Ecosystems: A Quantitative and Probabilistic Risk Assessment Perspective

Oil spills resulting from maritime accidents pose a poorly understood risk to the Arctic environment. We propose a novel probabilistic method to quantitatively assess these risks. Our method accounts for spatiotemporally varying population distributions, the spreading of oil, and seasonally varying species-specific exposure potential and sensitivity to oil. It quantifies risk with explicit uncertainty estimates, enables one to compare risks over large geographic areas, and produces information on a meaningful scale for decision-making. We demonstrate the method by assessing the short-term risks oil spills pose to polar bears, ringed seals, and walrus in the Kara Sea, the western part of the Northern Sea Route. The risks differ considerably between species, spatial locations, and seasons. Our results support current aspirations to ban heavy fuel oil in the Arctic but show that we should not underestimate the risks of lighter oils either, as these oils can pollute larger areas than heavier ones. Our results also highlight the importance of spatially explicit season-specific oil spill risk assessment in the Arctic and that environmental variability and the lack of data are a major source of uncertainty related to the oil spill impacts.


Extended materials and methods
In the following, we give a detailed description of the procedures used to estimate the parameters and functions included in the Kara Sea (KS) case study.

The OSRA calculations: solving the probability distribution for avgPPR
Here we describe in detail how the OSRA calculations summarized in the main paper and Figure 1 are done in practice. As an end product of OSRA, we want to solve the probability distribution of the expected oil spill impact conditional on a given route , season , oil type and volume of spilled oil which are the scenario variables that can be controlled by management ( Figure 1 in the main paper). This probability distribution can we written as Eq. S1 where we have used the conditional independence assumptions encoded in the DAG in Figure 1 of the main paper. If these conditional independence assumptions did not hold, the joint distribution of the OSRA factors could not be factorized as was done above, and we would need to solve their joint distribution directly. Note also that since PPR is defined to be a deterministic function of and , the probability distribution (PPR| , ) is a delta measure at ∫ ( , , ) ( , ) ∈ ( ) . Similarly is a deterministic function of and so that p( | , ) is a delta measure at ∫ (, ) ∈ ( , , , )̃. Moreover, the notation … in p( | , … ) denotes that the probability distribution for the accident probability should depend on factors not included in this current study. For this reason we did not estimate but fixed it to be a constant across the Kara Sea. By using the above conditional independence structure, we can solve p(OSI|d, t, o, V) by first solving each of the conditional probability distributions (Sections SI1.2-1.5) and after that marginalize over them. This can be done efficiently with Monte Carlo (MC) approximation as detailed in Algorithm 1 below.

Environmental covariates
The environmental covariates were chosen based on their importance to AMM habitat selection and oil spreading (1, 2; see Table S1). The raster maps of environmental covariates were derived from published scientific papers and open source data sets summarized in detail in Table S2. The environmental covariates, having different original resolution, were scaled to a lattice grid of 5 km. Temporal variation of sea ice concentration and sea surface salinity (SSS) were kept on original monthly resolution as bathymetry and distance to coast vary only spatially. The spatial and temporal resolutions were chosen based on the spatial and temporal accuracy of the data for species and environmental covariates. The full temporal extent of sea ice data was from 1996 to 2013, which includes the years when the species observations were recorded (Section SI1.3). The SSS data represents monthly averages in the study cells over the predicted values in years 1980-2000. Hence, it does not pose annual variation nor cover the latest years. However, the monthly averages capture the seasonal variation of SSS, and they describe well the area of fresh water and high nutrient impact in the KS (3).
In order to account for the uncertainty related to the seasonal variation in sea ice conditions in OSRA, we used their empirical distribution for years 2009-2013 to represent the annual and seasonal variability in the environment. The Arctic environment has been changing rapidly (4)(5)(6), which is why we restricted the data to include only the last five years from the sea ice data set. These years represent the present state of the Kara Sea ice conditions. SSS data are based on model predictions, which are made, among other covariates as well, according to river inflow from the large Siberian rivers, river Yenisei and river Ob. Due to a lack of inflow data, we could not cover the years 2009-2013 for SSS, and we needed to rely on monthly averages (6). We still think it is safe to use such temporally restricted SSS data in the species distribution models to represent the area of fresh water impact, where S3 high loads of nutrients increase the primary productivity in open sea season (7)(8)(9). Although, the average SSS over the Kara Sea may have changed since 2000, we assume that the area of fresh water impact and the seasonal pattern of how it forms in spring and weakens in autumn, have stayed similar to nowadays. Hence, in Algorithm 1, ( | ) corresponds to the empirical distribution of environmental covariates at season , and sampling of was performed as follows: 1. Sample a random month from all of the months corresponding to season during years 2009-2013. 2. Form a matrix of all environmental covariates in all raster cells in the Kara Sea in that random month. 3. Repeat steps 1 and 2 for all = 1, … , in Algorithm 1.

____________________________________________________________________________________
Note that the algorithm is generic in the sense that we can update the model for each individual conditional distribution without changing other components of the algorithm. This makes it easy to revise the analysis, for example, if more data are collected * ( ) collects the lattice grid cells through which route goes and | ( )| is the number of those cells. § As discussed in the main text, the accident probability is treated in this work as a constant throughout the Kara Sea. For this reason, we set , = 1/| ( )| so that avgPPR and avgOSI correspond to proportion of population at risk under the assumption that an accident happens as randomly selected location along the route. However, this assumption could be relaxed in the subsequent works so that , would be modeled as a random variable that depends on environmental conditions and possibly some other factors not included in the current model.

Species relative density
Species' relative densities were estimated with Mäkinen and Vanhatalo's species distribution model (SDM) (1), which links species count/occurrence observations to environmental conditions, spatial location and time. Here, we summarize their model and how it was used to produce the probabilistic species relative density predictions. For detailed description of the model, see (1).
Species distribution data. The data used to estimate species' relative densities consist of species observations and environmental data. The environmental data were the same as those described in Section SI1.2.
Species count and occurrence data are scarce in Arctic marine areas. Therefore, we collected data on polar bears, ringed seals and walrus in the KS from all available scientific articles and books (Table S3). These data consisted of species observations from the years 1996 to 2013. These observations were originally made from research vessels and ice breakers and reported in maps and tables summarizing when and where the surveys had taken place, and when, where and how many individuals of each species were observed. Some species observations were presence-only observations. The observations covered all seasons and thus, we were able to analyze the seasonally varying distribution of the species. The original publications did not report their study methods in detail. Hence, we did not know the probability of observing an individual or a group of individuals of a species, and had to account for this uncertainty in our model.
We digitized the maps and tables of species data in 5 km spatial resolution month-by-month. These resolutions were chosen based on the spatial and temporal accuracy of species and environmental covariate data. We compiled the species data with the environmental covariate data, resulting in a data table comprising a list of species observations (count, occurrence, absence) accompanied by information about spatial and temporal location and environmental covariates. The data contained in total 393 presence or positive valued count observations for polar bears, 231 for walrus and 1512 for ringed seals. Further, there were 8762, 9245 and 8446 records of absence, respectively.

Species distribution model (SDM).
The SDM followed the general Bayesian hierarchical spatial structure (10,11), which divides the model into an observation model, a process model and a parameter model. The observation model described the conditional distribution of the number of individuals, ( , , ), observed in a grid cell with coordinates (kilometers) and at time (months) during survey (in total three surveys of polar bears and walrus each, and one survey of ringed seals, Table S3) with a negative-Binomial distribution function, ( , , )| ( , , , ), ,~− ( ( , , , )+ , ).

Eq. S2
Here the latent function, ( , , , ), denotes the logarithm of the relative density of a species, , the vector of environmental covariates at grid cell at time , is a sampling effect of 'th survey and the overdispersion parameter accounting for spatially and temporally uncorrelated variation that is not explained by covariates, spatiotemporal random effect or sampling effect (to be described below). The spatiotemporally constant parameter for sampling effect adjusts for the variability in the abundance counts originating from the various sampling methodologies used by different surveys (that is, different data sources). It was modeled as independently and identically distributed Gaussian random variables, ~(0, 2 ), where 2 governs the variation in the sampling effects.
The log relative density corresponds to the process model and was modeled as ( , , , ) = + , + ( , ), Eq. S3 where is an intercept, = [ 1 , 2 … , ] is an ⨯ 1 vector of coefficients and ( , ) is a spatiotemporal random effect, which captures spatiotemporal variation that cannot be explained by the covariates. The vector of covariates, , , included all the environmental covariates and their squares, so the responses along environmental

S5
covariates were assumed to be quadratic. This is justifiable since conditions may be most favorable for the study species in the middle of the environmental gradients and thus, responses can follow a unimodal function. In addition to abiotic effects, we used the maximum a posteriori (MAP) estimate of relative density of ringed seal as a covariate in the polar bear model. We modeled the effect of ringed seals on polar bears with a Michaelis-Menten function, . It defines an asymptotic response between polar bear log relative density and ringed seal relative density, where denotes saturation and , half-saturation level. See (1) for justification and further details. Spatiotemporally varying random effect is given a Gaussian Process (GP) prior. GPs are a family of stochastic processes, which are flexible tools for modelling dependencies between observations in space, time and covariate space (10). We used a zero mean GP with a separable covariance function that is a product of squared exponential spatial and exponential temporal covariance functions where ST 2 is the process variance and , = 1, 2, 3 are the length-scale parameters governing how fast the correlation between ( , ) and ( ′, ′) decreases.
The last level of hierarchy is the parameter model, which defines the prior distributions for the parameters of the process model ( , , , covariance function parameters, etc.). The priors for parameters are summarized in (1).
We denote by Y= The posterior distribution was estimated with Markov chain Monte Carlo using the GPstuff toolbox (12) so that we sampled random samples of model parameters from the posterior distribution (see the work by Mäkinen and Vanhatalo (1) for further technical details). The convergence of Markov chains was analysed with the Gelman-Rubin Potential Scale Reduction Factor (PSRF). The models were validated with posterior predictive checks and cross-validation (13). As reported by Mäkinen  Next, the posterior samples and SDMs were used for the scenario predictions for the population distribution of the species in the KS. Since we did not have information on detection probability, the process level model does not model the expected abundance of a species as an absolute count of individuals, but as a relative density index. Hence, the latent function ( , , ) corresponds to the log relative density of a species in location s at time . However, although the available data does not estimate animal counts in an oiled area, it does estimate proportion of the KS population in that area. Since the relative density of a species is proportional to the absolute density, we can use the relative density to calculate the population distribution within the KS. Since the population sizes of polar bears, walrus and ringed seals show large seasonal variation in the KS (see Figures S1-S3), we use the S6 estimates of spring population sizes, which is typically the largest, to represent the actual sizes of the KS population. The population distribution is then Eq. S6 where tot is the area covering the KS and ̅̃ is a vector of average environmental covariates in spring. The proportion of the KS population under the oiled area, after a spill at location s, is then .

Eq. S7
In the implementation of OSRA, we need samples (see Algorithm 1) from its posterior predictive distribution Note, since each MC sample of , was calculated using the corresponding samples and , the samples of population distribution and the proportion of the population in an oiled area represent uncertainty in both environmental conditions ( , ) and knowledge of species distribution (parameters of SDM). This way the samples of ( , , , , ) account for the fact that, after marginalizing over the environmental conditions (, ) and ( , , , ) are not independent in their joint distribution. In practice, the above sampling steps can be done and samples stored in advance before running the whole Algorithm 1.

Oiled area
Oiled area at route points. Various physical and chemical processes such as spreading, evaporation, dispersion and dissolution determine the fate of oil in marine environment. These processes are affected both by the properties of spilled oil and by environmental conditions. In cold environments, also ice plays a critical role, as the oil can appear on ice, in ice, under ice, or trapped between blocks of ice (14,15). In recent years, there have been major advances in research and modeling related to spreading and fate of oil spills in ice conditions in the Arctic (see e.g. 14, [16][17][18][19][20]. However, the knowledge about oil spills in waters covered totally or partly by ice is still limited (16). Moreover, the current state-of-the-art oil spill fate and trajectory models typically require extensive data about bathymetry and weather and ice conditions, which are currently not available for many parts of the Arctic. Further, it is not usually feasible to conduct trajectory modeling with many potential spill locations, and therefore S7 the modeling is usually based only a few spill locations (e.g. 17,18,20). To overcome these issues, we applied the following approach.
Oil spreading in the sea is typically described by Fay's (21) equations, which have been modified over the years to predict areal coverage of oil in different environmental conditions. Here, we applied Fay-type equations developed by (2) for investigating the physical and chemical fate of oil in pack ice. For simplicity, we assumed that the release of oil is instantaneous and occurs in pack ice; and that the oil slick is circular throughout the duration of oil spreading. Moreover, we assumed that oil spreading occurs within a gravity-viscous regime (21), which dominates spreading in the longer term and is therefore the basis for most oil spill models. We also neglect the effect of evaporation because at cold water temperatures, it is predicted to have only a minor effect (see "Sensitivity analysis for the effect of evaporation" below).
With these assumptions, the basic Fay-type spreading model adjusted to different ice conditions describes the spill area, ( 2 ), as a function of time as where is the density of water ( 3 ), is the density of the oil ( 3 ), ∆ is the fractional buoyancy of oil ( − ), is the gravitational acceleration ( 2 ), is the volume of spilled oil ( 3 ), is the elapsed time ( ), is the is the viscosity correction factor and is the fraction of water covered by ice (%). When calculating the oiled areas, we fixed the model parameters to the values presented in Table S4. The slick sizes for light and extra heavy oils in different ice concentrations are shown in Figure S5.
Given the oil type and volume and after sampling from the distribution of environmental covariates (Section SI1.2), the construction of MC samples for , (see Algorithm 1) proceeded as follows: 1. Extract the ice concentration corresponding to 'th route cell from the matrix and use this ice concentration to calculate the the oiled area , . 2. Repeat step 1 for all ∈ ( ) as in Algorithm 1. 3. Repeat steps 1 and 2 for all = 1, … , in Algorithm 1.
The MC samples of environmental covariates used in calculating oiled areas were the same as those used to calculate population distribution (see Section SI1.3). This was then repeated for all oil types. In practice, the values for ( , , , ) at different ice conditions and oil types were first tabulated for different ice concentration values with 5% unit intervals to produce a look-up table for the above calculations. Since we fixed the model parameters, the MC samples only account for uncertainty due to environmental variability, but not for uncertainty in oil model parameters. Figure S4 demonstrates the overlap of a circular oil slick and species' densities in an example location along route 3.
Sensitivity analysis for the effect of evaporation. In order to estimate the effect of omitting evaporation from the estimates of oiled area, we calculated the theoretical maximum evaporation percentage of different oil types in the KS area using equations by (22). Since the evaporation of oil depends positively on water temperature (22), we assumed that the evaporation is highest during summer. Hence, we used the mean surface water temperature for summer months (July-September) in the KS area in the calculations (6°C). We selected Arabian light (API 31.3), Prudhoe Bay (28.5) and Bunker fuel (11.4) to represent light, medium, and heavy oil, respectively. It should S8 be noted that these same oil types were used in the spreading equation for consistency. The evaporation equations for these oil types were: where %EV is the percentage of oil evaporated, T is the sea surface temperature (°C) and t is time (min).
The evaporation percentages during summer months for light, medium and heavy oil after 14 days were 26.3, 26.0 and 5.7, respectively. If we assume a linear dependency between the volume of spilled oil and the area of the oil slick, we can see that although the inclusion of evaporation decreases the area of oil slick, especially for light and medium oils, the ratio between heavy and light or medium oil does not change substantially (Table S5). Hence, light oils have the potential to contaminate much larger areas than heavier ones, even if evaporation was taken into account, and therefore, our conclusions about the exposure potential of different oil types would not change significantly.

Exposure potential and sensitivity of species
We used expert elicitation and literature review (see Table S6) to acquire estimates of exposure potential ( , ) and sensitivity ( , ) for the AMM species included in the study (Figures S6-S8). These estimates and the uncertainty in them were encoded as probability distributions ( | , ) and ( | , ). For ringed seals, we used the estimates given by three independent experts in (23). In general, the experts assessed both the exposure potential and sensitivity of ringed seal to be relatively low under all accident scenarios (i.e. in all combinations of seasons and oil types). They estimated medium and heavy oils to be the most harmful oil types, and included very little variation in their estimates between seasons.
For polar bear and walrus, we conducted a literature review and used the elicitation tool built by Nevalainen et al. (23) to produce the estimates for these species. The probability distributions were initially developed by M. Nevalainen and later reviewed and revised in cooperation with the other authors. The literature search covered both scientific and gray literature. We used databases such as Scopus and Web of Science and keywords such as oil, oil spill, polar bear, walrus, habitat, behavior, toxicity, sensitivity, exposure etc. The search was limited to literature published in English.
With respect to sensitivity, walrus were assumed to be similar to ringed seals (24). Hence, the probability distribution for ringed seal sensitivity (23) was used also for walrus. The sensitivity of polar bears was assumed to be high with high certainty under all accident scenarios (based on, for example, (24) and (25)).
Both the exposure potential of polar bears as well as walrus were assigned a relatively wide probability distribution to indicate the lack of knowledge. The probability distribution for the exposure potential of walrus was assumed to resemble that of ringed seals (23,24), but more uncertainty was included in the estimates to recognize the inaccuracy related to such a generalization. Both ringed seals and walrus are skilled swimmers, which exploit ice and dive for food. However, walrus feeds on seafloor and therefore, may be more likely to come in contact with sunken oil on the seafloor. Hence, the exposure was highest to heavy oil, which was assumed to adhere both the seafloor and shoreline or ice edge, thus potentially smudging both the walrus' feeding and resting habitats. Nevertheless, walrus avoidance behavior related to oil has not been studied and therefore, the assigned probability distributions were wide under all accident scenarios. Polar bear was assigned an even wider exposure potential probability distribution than walrus to indicate the polar bear's highly dynamic lifestyle (26) and contradictory knowledge about their ability to avoid oil in their habitat

S9
(see (27) vs. (28)). The polar bear's exposure potential was highest for heavy and medium oils, as these oils can smother ice and other habitats and accumulate in openings in ice, the polar bear's hunting habitat. Due to lack of knowledge that would suggest there to be seasonal differences in exposure potential, the same probability distributions were used for all seasons.
As discussed by Nevalainen et al. (23) and O'Hagan (29), assessing the goodness of the elicited probability distributions for sensitivity and exposure potential is hard and quantitative assessment is virtually impossible since there are no data to inform about the true sensitivity and exposure potential. As our interviewees were (natural) scientists, it is likely they were comfortable working with probability distributions. Moreover, the interactive graphics in the elicitation tool was intended to help the experts to assess whether their beliefs were accurately presented by the probability distributions they assigned. Such a feature has been found useful in many expert elicitation studies so we believe that the results of the elicitation reflect our experts' knowledge with sufficient accuracy despite the shortcomings discussed here. Since we do not have empirical data on exposure potential and sensitivity, it is impossible to quantitatively measure how well the results reflect reality. Hence, we can only compare the elicited results to prior knowledge about oil spill impacts to assess the goodness of the probability distributions. The expert assessments in this study are in line with literature and this enhances their credibility. For more in depth discussion on the potential biases and errors in the sensitivity and exposure potential see the work by Nevalainen et al. (23).
Given the oil type and season the MC samples of sensitivity , and exposure potential , , to be used in Algorithm 1, were created by sampling respectively from ( | , ) and ( | , ). S10 Figure S1. The expected relative densities of polar bears in the Kara Sea with a monthly resolution in 2009-2013. The inter and intra annual variation in expected relative density arises from variation in environmental covariates. S11 Figure S2. The expected relative densities of walrus in the Kara Sea in a monthly resolution in 2009-2013. The inter and intra annual variation in expected relative density arises from variation in environmental covariates. S12 Figure S3. The expected relative densities of ringed seals in the Kara Sea in a monthly resolution in 2009-2013. The inter and intra annual variation in expected relative density arises from variation in environmental covariates. S13 Figure S4. Example of the overlap of the oil slick (black circle) and the species densities at high (left column) and low (right column) ice concentrations. Upper three rows: polar bear, walrus, and ringed seal, respectively. Bottom row: Ice concentration in the KS. The red square shows the location used as an example in the figures in upper rows; the blue line represents a shipping route along the KS (route 3).   Figure S5. Oil spreading under different ice cover conditions within 14 days after an oil spill. A. Light oil. B. Extra heavy oil. Table S4. Parameter values used in oil spreading calculations. Note that for oil parameters, we used a fixed value representing a typical oil in a given class (shown in bold). The range of typical values for each class is shown in parentheses.

Abbreviation Explanation avgPPR
Expected proportion of the population (of a given species) that occurs in the oiled area (i.e. the Proportion of Population at Risk), when an accident occurs at a random location along a route PPR Route-specific risk scaled by route length avgPPR_cmb Combined avgPPR, defined as the sum of individual avgPPRs of the three AMM species avgOSI Expected Oil Spill Impact, i.e. the proportion of the population that dies due to oiling, when an accident occurs at a random location along a route, defined as avgOSI=avgPPR×BI OSI Route-specific oil spill impact scaled by route length avgOSI_cmb Combined avgOSI, defined as the sum of individual avgOSIs of the three AMM species BI Biological Impact (of a given species), defined as BI=exposure potential (of a given species) × sensitivity (of a given species) S22 Figure S9. AvgPPR for AMM species along routes R1-R5 with different oil types (rows) and seasons. Note the different scales (but same units) on the y-axes. Circle denotes median, thick line the 25% to 75% quantile, and thin line the 5% and 95% quantile.
S23 Figure S10. PPR for AMM species along routes R1-R5 with different oil types (rows) and seasons. Note the different scales (but same units) on the y-axes. Circle denotes median, thick line the 25% to 75% quantile, and thin line the 5% and 95% quantile.
S24 Figure S11. AvgOSI for AMM species along routes R1-R5 with different oil types (rows) and seasons. Note the different scales (but same units) on the y-axes. Circle denotes median, thick line the 25% to 75% quantile, and thin line the 5% and 95% quantile.
S25 Figure S12. OSI for AMM species along routes R1-R5 with different oil types (rows) and seasons. Note the different scales (but same units) on the y-axes. Circle denotes median, thick line the 25% to 75% quantile, and thin line the 5% and 95% quantile.