Species and habitats in danger: estimating the relative risk posed by oil spills in the northern Baltic Sea

Largescale oil spills can have adverse eff ects on biodiversity in coastal areas where maritime oil transportation is intense. In this article we conducted a spatial risk assessment to study the risk that potential tanker accidents pose to threatened habitat types and species living in the northern Baltic Sea, which has witnessed a rapid increase in maritime oil transportation within the past two decades. We applied a probabilistic method, which combines three components: a Bayesian network describing tanker accidents and uncertainties related to them, probabilistic maps showing the movement of oil, and a database of threatened species and habitats in the area. The results suggest that spatial risk posed by oil spills varies across the area, and does not correspond, for example, to the frequency of accidents in a given area. The relative risk is highest for seashore meadows, which is important to take into account when managing these habitats. Our analysis underlines the importance of a thorough risk assessment, which is not only based solely on one or two specifi c factors such as accident probabilities or the trajectories of spilled oil but also contains as broad a view of the consequences as possible. We believe that the probabilistic methodology applied in the study will be of high interest to people who have to cope with uncertainties typical for environmental risk assessment and management.


INTRODUCTION
Various human-induced stressors like habitat change, overexploitation, invasive species, and nutrient enrichment have led to a signifi cant loss of biodiversity in marine and coastal ecosystems across the globe (Millennium Ecosystem Assessment 2005 ). Another severe environmental problem aff ecting these ecosystems is accidental or chronic oil pollution (NRC 2003 ). Although the total number of oil spills resulting from tanker accidents has been declining globally in recent decades (Musk 2012 ), accidents like Exxon Valdez in 1989, Erika in 1999, and Prestige in 2002 have had major negative impacts on the environment (Peterson et al. 2003, Cadiou et al. 2004, Penela-Arenaz et al. 2009 ).
Oil can alter marine and coastal ecosystems in various ways. The negative eff ects on organisms usually arise from physical smothering, toxicity of aromatic hydrocarbons, habitat modifi cation, or the combination of these (Albers 2003  Exposure to oil can lead to acute death as well as varying kinds of sublethal eff ects that decrease the fi tness of an individual, and even the structure and function of whole communities may change due to such exposure (NRC 2003 ). The effects of oil on the biota have been studied extensively in laboratory se ings as well as a er major oil accidents (for examples see reviews by Paine et al. 1996, Albers 2003, NRC 2003, Peterson et al. 2003, Penela-Arenaz et al. 2009 ), and sophisticated oil spill impact models have been developed (e.g., SIMAP;French McCay 2003. However, due to the irregular nature of oil spills and the typical variability in oil spill impacts it is diffi cult to arrive at any straightforward conclusions related to the eff ects of oil spills and thus the risk posed by oil spills (NRC 2003 ). Ecological, or environmental, risk assessment (ERA) aims at predicting the impacts of various stressors on the environment. ERA has two interlacing core concepts: decision-making and probability (Suter 2007 ), as the main purpose of ERA is to support decision-making under uncertainty (Burgman 2005 ). Although ERA was originally developed for assessing the eff ects of specifi c hazardous chemicals on organisms living in specifi c areas, it nowadays covers various types of stressors and ecological consequences (Suter 2007 , Ayre andLandis 2012 ). Also the term "risk" can be defi ned and quantifi ed in multiple ways (e.g., Kaplan andGarrick 1981 , Fischhoff et al. 1984 ), but in oil spill-related research a suitable choice is the widely used defi nition that risk is a combination of the probability of an adverse event and the consequences of the event (Burgman 2005, International Organization for Standardization 2009. The importance of spatial risk assessment has been recognized in an oil spill management context (Frazão Santos et al. 2013 b ). For instance, spatial risk assessment can help with the deployment of oil combating equipment along shorelines and suggest areas where special a ention to oil combating contingency planning is needed. In addition, the results of oil spill risk assessment can be used to guide marine spatial planning (Frazão Santos et al. 2013 b ), including, for example, the planning of environmentally sound shipping routes (Soomere et al. 2011 a ).
In this article we apply a probabilistic method to assess the spatial risk that oil tanker accidents pose to coastal ecosystems. Our approach encompasses three spatial components that are relevant to a comprehensive risk assessment: (1) Varying probabilities and locations of accidents, (2) The stochastic spreading of oil, and (3) The coastal values at stake. Our work is based on the methodology for a technical procedure to combine the three components developed by Jolma et al. ( 2014 ).
We use the northern Baltic Sea for a case study to demonstrate our approach. Despite the high numbers of oil tankers and other vessels operating in the area (see HELCOM 2010 ), there are only a few previous studies of oil transportation risks with at least some kind of spatiality. Kujala et al. ( 2009), Goerlandt et al. ( 2012, Hänninen et al. ( 2012 ) and Lehikoinen et al. ( 2015 ) studied accident frequencies and subsequent oil spills in some specifi c locations, and Soomere et al. ( 2011 a , b ) applied a Lagrangian-type transport model to fi nd an environmentally optimal fairway, which would minimize pollution on shorelines in the area. Further, Helle et al. ( 2011 ) analyzed the threat posed by oil spills in the Western Gulf of Finland by combining the spreading of oil with species occurrences, but their analysis was limited to one spill location and a few species. Hence, comprehensive analyses that consider all three spatial aspects relevant to oil spill risk assessment are still lacking.
We concentrate on threatened species and habitats, as they are usually already strongly aff ected by other human activities, and, as argued by Ihaksi et al. ( 2011 ), the recovery of threatened and rare species is slow and uncertain compared to common species. Further, there are relatively comprehensive databases on the occurrences of threatened species and habitats in Finland, which enables us to test and apply the data in a spatial risk assessment context. First we analyzed the overall spatial risk in order to identify areas that are most prone to oil, a er which we conducted a more detailed analysis to fi nd out which habitat types and species in the study area are most vulnerable to oil, that is, have the highest relative oil spill risk scores.

Study area: the Gulf of Finland and the Finnish Archipelago Sea
The Gulf of Finland (herea er GoF) and the Finnish Archipelago Sea (AS) are located in www.esajournals.org HELLE ET AL. the northern Baltic Sea, one of the largest brackish water bodies in the world (Fig. 1 ).
The GoF is a straight continuation of the Baltic Proper, and it borders Finland, Estonia, and Russia. It is a narrow and shallow gulf 44-135km wide with an average depth of 37 m (Alenius et al. 1998 ). The coastline is indented with numerous islands and islets especially in the northern side of the gulf. The salinity of its waters vary from 0 PSU (practical salinity unit) in the east to 7 PSU in the west, the gradient being formed as a balance between more saline water protruding from the Baltic Proper and freshwater fl owing from rivers (Alenius et al. 1998 ). The Finnish AS is located westwards from the GoF, between the Baltic Proper, the Åland Sea, and the Gulf of Bothnia. The defi ning feature of the area is the vast archipelago, which is comprised of over 22,000 islands (Granö et al. 1999 ). The mean water depth of the AS is 23 m, and the mean salinity is 6-7 PSU (Haapala and Alenius 1994 ). In winter, the GoF as well as the AS freeze usually at least partly, the average ice period lasts between 70-110 and 40-130 d, respectively (Seinä and Peltola 1991 ).
The brackish water together with regular icecover makes the study area an exceptional place among the world ' s marine ecosystems. As in the Baltic Sea in general, the aquatic biota is a mixture of marine and freshwater species (Elmgren and Hill 1997 ), and the coastline is a mosaic of varying biotopes from steep cliff s to sandy beaches and sheltered bays covered with vegetation. Due to the last ice age, the GoF and AS have features that are not present anywhere else on the Finnish coast. For example, the Hankoniemi area in the western GoF and the island of Jurmo in the AS are continuations of the Salpausselkä ridge systems, and contain biotopes atypical for the area. In addition, due to the postglacial upli new habitats like glo-lakes and fl ada-lakes are still forming. Especially in the AS area the archipelago further diversifi es the environment, and the area is characterized by high geodiversity, which together with traditional forms of land use (e.g., pasturage) have led to high biodiversity of coastal areas (Granö et al. 1999 ).
In recent decades, maritime oil transportation has increased substantially in the study area. The growth has been distinct especially in the GoF, where the estimated volume of transported oil has increased eightfold since 1995 being approximately 160 million tons in 2013 (Finnish Environment Institute 2013 ). The main cause for this development is increased Russian oil exports, reinforced by the construction of new oil terminals and the modernization of the old ones (Knudsen 2010 ). There are altogether 15 oil terminals in the GoF, and the biggest oil terminal in the Baltic Sea is located in Primorsk (Koivisto) in the eastern part of the area (Brunila and Storgård 2012 ). The current oil terminal capacity in the GoF is almost 250 million tons. Concomitantly, the AS harbors one oil terminal (Naantali).
Intensifi ed oil transportation has raised concerns about the risks of major oil accidents in the study area. A realistic worst-case scenario for both areas is estimated to be a collision between a vessel and an oil tanker carrying crude oil. Taking typical tanker sizes into account, this could result in a spill of 30,000 tons and 15,000 tons in the GoF and the AS, respectively (Hietala and Lampela 2007 ).

Risk assessment framework
The risk assessment is based on an integration of three elements (Fig. 2 ), which all have a spatial dimension: (1) A Bayesian network model which defi nes the characteristics of possible tanker accidents in the study area in a probabilistic manner, (2) Probabilistic maps describing the movement of oil a er an accident, and (3) Data on the locations of threatened species and habitats along the Finnish shoreline. The methodology to combine these three elements was originally demonstrated by Jolma et al. ( 2014 ).
Bayesian network: uncertainty related to the accidents.-The fi rst element is a Bayesian network model describing the important features of oil accidents and the uncertainty related to them. Bayesian networks (BNs) are graphical models that represent probabilistic relationships in a set of variables relevant to a problem in a causal manner. The variables (i.e., nodes) together with the edges (i.e., arrows) between them form a directed acyclic graph (DAG). A variable that is dependent on other variables is called a child, whereas a variable on which other variable is dependent is called a parent (Jensen and Nielsen 2007 ). Although some BN so ware can also handle continuous (Gaussian) distributions, a standard procedure is the discretization of the variables. The dependency between the variables is described with a conditional probability table (CPT), which contains the information of the probabilities of a child being in a certain state given the states of its parents. By applying the rules of probability calculus and Bayes' theorem, BNs can calculate not only from causes to consequences but also vice versa, that is, from consequences to causes. Thus, when information is fed into any of the variables in the model (e.g., if a variable is locked Fig. 2 . The components of the risk assessment. Ovals represent probabilistic variables, and variables inside black dashed line are included in the BN used in the study. The variables shown in gray are input for oil spill model Spillmod, which is used to produce probabilistic maps describing the movement of oil. The probability distributions of the input variables are further used to weight the maps. The fi nal calculation combines the information of probabilistic maps with the species and habitat data.
to a certain state), the probability distributions of other variables change accordingly depending on the CPTs.
In addition to other fi elds of science, BNs are increasingly applied in environmental science and ecology (see, for example, reviews by Mc-Cann et al. ( 2006 ) and Aguilera et al. ( 2011 ) ). Bayesian networks off er a convenient way for reasoning under uncertainty, and they have many features that make them useful tools also in risk assessment, as they, for example, model causality explicitly and can combine diff erent types of evidence such as subjective expert knowledge and statistical data (Fenton and Neil 2013 ). In recent years, BNs have also been applied to spatial risk assessments, for example, by Stelzenmüller et al. ( 2010 ) and Ayre and Landis ( 2012 ).
We used a BN to describe diff erent tanker accident scenarios in a probabilistic manner. The BN applied contains fi ve variables (Fig. 2 , Table 1 ). In addition, a sixth variable ("Season") is included in the risk calculations as a probabilistic variable, although it is not linked to any other variable in the BN. The variables in the BN represent the main factors relevant to the accident, and the focal assumption is that a spill has occurred.
The variable "Spill size" describes the size of the spill in tons. It is dependent on "Tanker size," "Accident type," and "Oil type". The fi rst one describes the size of the tanker involved in the accident, and the second one the type of the accident, that is, grounding or collision. The third one represents the type of spilled oil, and it has three states: light, medium, and heavy. The variables are dependent on "Location," which represents the location of the spill. Eight potential spill locations (fi ve in the GoF and three in the AS), distributed relatively evenly across the study area, were included in the analysis (Fig. 1 ). The locations were chosen to match the crossing or merging points of important maritime routes in the area. The variables and their discretized states are described in Table 1 .
The prior distributions for "Location," "Accident type," and "Season" were gained from the HELCOM (Helsinki Commission) statistics of shipping accidents in the Baltic Sea (HELCOM 2014 ). In order to compare the number of accidents between diff erent locations within the study area, a circle with a radius of 20 km was drawn around each spill location (see Fig. 1 ); the size of the circle was chosen so that the circles did not overlap. Within the circles, all groundings and collisions from the HELCOM database from 1990 onwards were counted and divided according to the type of the accident and season. However, as there have been only a few tanker accidents in the study area, we included all vessel types in the analysis, that is, we assume that oil tankers behave in the same way as other vessel types reported in the database (mainly cargo and passenger vessels in addition to tankers). We also took into account the ratio of tankers to other vessels operating in the study area.
Similarly to the number of accidents, the probability distribution for "Season" was formed by using all available data, that is, by assuming that there are not any diff erences between areas due to the season.
The probability distributions for "Tanker size" and "Oil type" were obtained from Portnet statistics ( www.portnet.fi , year 2010). For these variables the spatial resolution was fairly coarse, that is, the GoF was treated as one unity and the AS as another one. For the spill location, AS3 the probability distribution was calculated as a weighted average of the GoF and AS distributions. The weighting was based on the volume of transported oil (oil type) or the number of vessels visiting the area (tanker size).
"Spill size" was calculated with lognormal distribution with the parameters mean (m) and standard deviation (SD). The mean outfl ow val- ues for both accident types were obtained by applying the method described by Montewka et al. ( 2010 ). More specifi cally, the mean outfl ow given there is a spill was calculated by using the outfl ow values for bo om and side damages reported by Herbert Engineering Corp. ( 1998 ), a er which a regression analysis was conducted with calculated outfl ow values against the reported tanker sizes. We converted the values reported in cubic meters to tons by using the specifi c gravities of 0.823 (light oil), 0.834 (medium oil), and 0.933 (heavy oil) kg/ dm 3 . The mean outfl ow values for each tanker size ("Tanker dwt") were used as m in the lognormal distribution, whereas SD was iterated by se ing the probability for unrealistic large spills (i.e., spills larger than the size of the tanker) to 0.0005. Probabilistic maps: spreading of oil.-The second component of the analysis is the probabilistic spatial data showing the movement of oil a er a tanker accident has occurred. The study area was covered by a grid, the cells of which had dimensions of 3 × 3 km, and the data describe the probability that there is oil in a certain grid cell within 240 h a er an accident given the characteristics of the leak and the prevailing weather conditions during and a er the spill. The data were produced by the Finnish Environment Institute with the oil spill model Spillmod (Ovsienko 2002 ). The main idea of the modeling was to compute the spreading of oil spill (with a duration of 8 h) taking place in a given location in the fi rst hour of a given month, in the second hour of a given month etc., and following the fate of the oil slick 240 h a er the release. The model output is a spatial data set, which describes the "hits" of the oil slick in a specifi c grid cell relative to all possible trajectories originating from the given accident location.
A separate data set was produced for each of the combinations of the relevant variables, which included location (eight alternatives, see Table 1 ), oil type (three), spill size (six), and season (three). The weather data used in the analysis comprise data from the beginning of March to the end of November for 1996-2001. Thus, there were altogether 2592 data sets (each consisting of 2184 or 2208 separate trajectory simulations) describing all the possible combinations of the four variables related to the spill. The years were combined to describe the average situation in a given season. A more detailed description of the data can be found in Jolma et al. ( 2014 ).
Habitat and species data.-The third component of the risk assessment that we carried out is the data of the occurrences of threatened habitat types and species living along the coastline. The data were gathered from diff erent sources. For threatened habitat types, the most important ones were the INSPIRE1 and MH SutiGis databases (maintained by the Finnish Environment Institute and Metsähallitus Parks and Wildlife Finland, respectively), whereas for threatened species the key sources were the Her a and Tiira databases (maintained by the Finnish Environmental Administration and Birdlife Finland, respectively). The data were stored in a separate database, which included the names and locations of habitats and species (polygons for habitats, points and polygons for species), their conservation status, and the area of certain habitat patch (habitats) or the number of individuals in or the area of a specifi c occurrence (species). The study area was delineated according to the national typology of coastal waters (Kangas et al. 2003 ). In the GoF all Finnish waters were included in the analysis, whereas the AS area was defi ned to include Finnish waters without the Åland Islands.
Several criteria were used to select species and habitats for the analysis. For habitats, the primary selection was based on the conservation status and the potential for the habitat type to become oiled a er a spill. Thus, each habitat type in the analysis (1) had the status of "Critically endangered" (CR), "Endangered" (EN), or "Vulnerable" (VU) according to the assessment of threatened habitat types in Finland (Raunio et al. 2008 ), and (2) was considered to be threatened by oil by a panel of experts. The experts were selected based on their knowledge of coastal ecology in the study area. The panel included mainly fi eld ecologists and conservation biologists, who were able to provide estimates, for example, of habitat representativeness and species behavior. Further, each separate habitat patch in the analysis was located within 100 m from the shoreline and was estimated to exhibit typical features of a certain habitat type by the www.esajournals.org HELLE ET AL.
experts. By applying these criteria, 8 habitat types and 610 separate habitat patches were included in the analysis (see Appendix S1 ). Underwater habitats were excluded because we have fragmented and incomplete knowledge of their conditions.
Each species included in the analysis had the following characteristics: (1) The IUCN classifi cation status of the species was CR, EN, or VU according to the most recent assessment of threatened species in Finland (Rassi et al. 2010 ), and (2) the reported primary habitat of the species was either the Baltic Sea or the shores of the Baltic Sea (Rassi et al. 2010 ). Similar to habitats, only occurrences within 100 m from the shoreline were taken into account. Altogether 75 species and 595 occurrences were included in the analysis, the majority of the species being bu erfl ies or moths (34.6%), vascular plants (17.3%), and beetles (17.3%) ( Appendix S2 ).
Also other a ributes relevant to the analysis were stored in the database. The most important is the VAL number, a numerical value given to a specifi c IUCN class (Ihaksi et al. 2011 ). The VAL describes the value of each IUCN class compared to other IUCN classes, and it also takes into account whether the species or the habitat is either a directive species or habitat (mentioned in annexes of the European Union ' s Bird Directive or the Habitat Directive) or is a so-called national responsibility species or habitat of Finland. The la er include species and habitats for which a large proportion of their total population or occurrences is located in Finland, and for which Finland can thus be seen to have an international responsibility (Rassi et al. 2001, Raunio et al. 2008 , for the concept see Schmeller et al. 2008 ). The relative values for each class were given by a panel of experts, who were chosen based on their knowledge of conservation legislation and management related to national and European species and habitats. The lowest value (0.2) was assigned to VU species and the highest (1) for species and habitats that were classifi ed as CR and were both directive and responsibility species/habitats (Table 2 ). Other relevant information for a species was the mortality value (MOR), which describes a species' sensitivity to oil ( Appendix S2 ). The MOR is defi ned as the proportion of a population that would die due to oil exposure if all individuals were in contact with oil (see Ihaksi et al. 2011 ). For habitats, a similar index is recovery potential value (REPH), which describes the probability that a habitat would recover within 5 yr a er oil exposure to a preexposure state defi ned by Raunio et al. ( 2008 ), if cleaned up properly (R. Venesjärvi, I. Helle and A. Jolma, unpublished manuscript; Appendix S1 ). More detailed information on survey and VAL, MOR, and REPH values can be found in Ihaksi et al. ( 2011 ) and Venesjärvi et al. (unpublished manuscript). The data have also been used in a web-based oil spill response tool designed to guide operational decision-making regarding shoreline combating operations (Altartouri et al. 2013 ).

Applying the risk methodology
When applying the defi nition that risk is the product of the probability of the incident and the subsequent consequences, the risk calculation procedure includes two steps. The fi rst step is to calculate the probability that a certain grid cell would become oiled a er an oil spill. This is done by weighting the oil spill probability maps by the probabilities describing the characteristic of a spill calculated in the BN (see Jolma et al. 2014 ). The second step is to calculate the value of each grid cell in order to be able to estimate the potential loss related to oiling (i.e., the consequence of the incident), a er which the fi nal risk score can be calculated.
Spatial risk.-We applied the basic concept that the risk score of the cell i is the product of the cell ' s probability of oiling (Pr(oil) i ) and the value of the cell (Value i ), that is, In order to compare the risk when, for example, a spill is supposed to take place in a certain location, the total oil spill risk score was calculated as the sum of risk scores of all cells, that is, where n is the number of cells aff ected by oil, and Pr(oil) i and Value i are the same as in Eq. 1 .
The value of each cell was calculated as the sum of the conservation values (VALs) of habitats and species present in the cell (Table 2 ). However, as some occurrences of some species had multiple observations, each species was calculated only once per cell.
Habitat-specifi c risk.-As we also wanted to compare the risk level that diff erent habitat types experience, we conducted a habitat-specifi c analysis, which takes into account both the probability of exposure and the area of exposed habitat, the la er expressing the importance of a certain habitat patch to the habitat type as a whole. Thus, the habitat-specifi c risk was calculated as (3) where n is the number of cells in which habitat H is present, Pr(oil) i is the probability that cell i becomes oiled, and propH i is the proportion of H present in cell i (i.e., the area of H in cell i compared to the total area of H in the study area).
However, as the malignancy of oil exposure varies between the habitats (Jensen et al. 1998 ), we also calculated a risk value that takes this into account. A more comprehensive habitat risk was then calculated as (4) where n , Pr(oil) i and propH i are the same as in Eq. 3 , and REPH H describes recovery potential for the habitat H .
In order to make the results more easily interpretable, all risk scores were scaled to the interval 0-1.
Species-specifi c risk.-The method to calculate the species-specifi c risk was similar to the habitat-specifi c risk calculation with some modifi cations. The optimal approach would have been to use the number of individuals or the area covered by the species to express the size of an occurrence. However, the data included many observations, where only one or two individuals were recorded, although it is realistic to assume that there are at least some more individuals present in the population. We, thus, applied two approaches: for easily observed species (birds and the majority of plants) we used the absolute numbers or the area reported in the database, whereas for less visible species we applied an adjusted fi gure. For these species, all observations of 10 individuals max were given numerical value of 10, whereas for observations greater than this absolute numbers were applied.
The species-specifi c risk was thus calculated as (5) where n is the number of cells aff ected by oil, in which species S is present, Pr(oil) i is the probability that cell i would become oiled, and propS i is the value of cell i , calculated as the proportion of species S present in cell i .
Similar to habitats, we also calculated a more comprehensive risk score, which takes into account also the mortality of individuals due to oil, not just the probability of exposure. The risk was calculated as (6) where n , Pr(oil) i and propS i are same as in Eq. 5 , and MOR S is the mortality value of the species S .
In order to make the results more easily interpretable, all risk scores were scaled to the interval 0-1.

Spatial risk
The probability distributions of the variables in the basic state of the BN are shown in Table 3 . The probability of oil spills is highest in the eastern GoF (G5 and G4), but the central GoF (G2) has also a reasonable likelihood to witness a spill (see also Fig. 4). Smaller spills seem to be much more probable than larger spills.
To some extent, ecological values cover the whole study area, but certain hotspots can be recognized, especially in the western GoF and in the southern AS (Fig. 3 a). When the BN is in the basic state and only the spreading of oil is considered, the highest probabilities can be found in the locations of accidents in the eastern GoF (Fig. 3 b). The fi gure also shows that in the study area oil typically spreads to the north-eastern direction. Furthermore, it seems that in AS the dense archipelago blocks oil from spreading to the inner parts of the area. When the two components are combined in the fi nal risk assessment, that is, each cell has a risk score calculated with Eq. 1 , risk seems to be highest in the outer archipelago of the eastern GoF, in the Hankoniemi area in the western GoF and in the southern part of the AS (Fig. 3 c).
We can test specifi c scenarios by instantiating one or more variables in the BN. We present a case, where the location of the spill is known, but all other variables are still uncertain. Fig. 4 shows how the relative risk (i.e., risk scores are scaled in relation to the highest value) varies depending on what components of risk we are interested in. As shown already in Table 3 , accident frequencies are highest in the eastern part of the study area (G4 and G5). However, when we assume that a spill happens in a certain location, and we therefore concentrate only on the spreading of oil and the ecological values for being oiled, the highest relative total risk scores RiskTOT (Eq. 2 ) are associated with locations AS1 and AS2 in the Archipelago Sea. This seems to be explained by the high number of ecological values in the area (Fig. 3 a). In the GoF, the highest risk scores are associated with locations G1 and G3, which is probably explained by both ecological values in the vicinity and the potential for oil to pollute large areas. The lowest values are related to accidents occurring in the easternmost GoF (G5), that is, the area, which has the highest accident frequencies.
Finally, when we multiply the location-specifi c RiskTOT scores with the probabilities of spills to happen in certain locations, we can compare different accident locations in a way that takes into account all three spatial components of oil spills (Fig. 4 ). This inspection shows highest risks for areas AS2 and G4, the former having high polluting potential and the la er having relatively high accident frequency.

Habitat-specifi c risk
The extent of risk posed by oil spills is different for the various habitats (Fig. 5 ). If only the exposure to oil is taken into account (RiskH), coastal sand beaches and seashore meadows seem to be especially at risk. When also the recovery potential of habitats is acknowledged (RiskCH), seashore meadows and Alnus glutinosa swamps, which have low recovery potential,

Species-specifi c risk
There does not seem to be any clear pa ern when species-specifi c risk scores are considered, and, for example, species with the highest Fig. 3 . Spatial distribution of ecological values (a); Spreading of oil when the BN is in the basic state (b); Spatial risk when ecological values are combined with the probability of oiling (c). Increasing values in a and c are represented with colors light green < light blue < purple < red < black. Note that in b, the contour lines crossing land area in some places is an artifact due to original data, which is a square grid of 3 × 3 km cells. relative risk include species from various taxonomic groups (Table 4 ). The primary habitat for the majority of these species is sandy beaches, which is logical as these habitats rank highest also in habitat-specifi c analysis (Fig. 5 ). RiskS and RiskCS seem to encompass almost similar information, which results from the fact that the estimated MOR values vary only a li le between species ( Appendix S2 ). For example, from 20 species with the highest RiskS scores 18 can be found also in the top 20 list of the highest RiskCS scores (Table 4 ). Some species like fi eld restharrow ( Ononis arvensis ) and sand couch ( Elymus farctus subsp. boreoatlanticus ) have high relative exposure scores (RiskS) but due to their somewhat lower MOR values, their comprehensive risk scores (RiskCS) are lower. Risk scores for all species included in the analysis are presented in Appendix S2 .

DISCUSSION
Oil spill risks exhibit a clear spatiality in the Gulf of Finland and the Archipelago Sea. The  Fig. 1 ) when the components of risk are assessed separately or in a combined manner. Black columns: Only accident frequencies are considered; White columns: It is assumed that an accident happens in a certain location, and only the dri ing of oil and ecological values at risk are taken into account; Striped columns: Accident frequencies (A) are combined with dri ing of oil (D) and ecological values (V). Note that as there have not been any reported accidents in AS 3, relative accident frequency as well as combined risk is 0 for this location. Fig. 5 . Relative risk scores of threatened habitat types, the calculation is based on the exposure to oil (RiskH, black columns), or the exposure and the recovery potential (Risk CH , striped columns). Habitat types: 1-Coastal sand beaches; 2-Seashore meadows; 3-Coastal dynes; 4-Alnus glutinosa swamps; 5-Calcareous rock outcrops on seashores; 6-Glo-lakes (coastal lagoons); 7-Flada-lakes (coastal lagoons); 8-Herb-rich forests with broadleaved deciduous trees.
results suggest that there are distinct diff erences in the overall spatial risk as well as in risks between individual habitats and species.
The spatial risk seems to be highest in the outer archipelago and in certain hotspot areas along the coastline. It is interesting that although the probability of tanker accidents is the highest in the eastern part of the study area, the situation is diff erent when the behavior of oil and the locations of threatened species and habitats are taken into account. If we consider only the pollution potential, that is, the stochastic dri ing of oil together with the ecological values exposed by oil, the Archipelago Sea seems to be the riskiest area. Accidents in the eastern part of the study area have the lowest pollution potential. Although this fi nding is at least partly explained by the fact that only ecological data from the Finnish coast were included in the analysis, it also refl ects the uneven distribution of these values. Furthermore, the result is in accordance with the work of Delpeche-Ellmann and Soomere ( 2013 ), who studied, how diff erent Marine Protected Areas (MPAs) in the GoF would be hit by currentdriven pollution originating from hypothetical sources along the main fairway. Their analysis showed that the length of fairways having a potential to pollute a MPA located in eastern GoF covers almost the whole GoF, whereas MPAs in the central and western GoF would be hit by the pollution originating mainly from the same areas, not from the eastern GoF. Oil accidents in the east thus seem to be more local than in other areas in the GoF.
Habitat-specifi c analysis suggests that the relative risk is high especially for seashore meadows, as they have high exposure probability (RiskH) but low recovery potential also resulting in a high comprehensive risk value (RiskCH). This is an alarming fi nding, as seashore meadows are already suff ering heavily from other anthropogenic stressors such as overgrowing, eutrophication and the construction of waterways, and are already classifi ed as critically endangered (Raunio et al. 2008 ). Alnus glutinosa swamps have not only somewhat lower exposure probabilities but also very low recovery potential. Their safeguarding from oil is thus essential, as even with careful clean-up the prespill status of the ecosystem may be impossible to achieve. Sandy beaches and coastal dynes, in turn, have high relative exposure probabilities, but their recovery potential a er appropriate clean-up is also fairly high. Yet, it is important to notice that they also harbor many threatened species living in close contact with sand, that is, organisms that may be harmed by clean-up operations such as the scraping of the uppermost layers of the sediment or pressurized washing methods (Paine et al. 1996, Pezeshki et al. 2000, Peterson et al. 2003, which makes clean-up a challenging task.
The importance of seashore meadows and sandy beaches are also highlighted when speciesspecifi c risks are considered, as the majority of the species that exhibit highest relative risk values live in these habitats.
Although coastal lagoons like glo-lakes and fl ada-lakes have fairly high RiskCH values, in reality they may be partly or even entirely safe from oil, as their connection to the sea is either weak or absent. This holds true also for some other habitats like coastal forests and dynes, which usually can be found in the upper parts of shores. However, they were included in the analysis, as they become exposed to oil in stormy weather or when the sea level is high. Furthermore, all these habitats have fairly low recovery potentials, and protective measures applied before any oil has reached the shoreline have an important role in mitigating the risk.
As pointed out by Frazão Santos et al. ( 2013 b ), spatial risk assessment can off er valuable information for the development of oil spill management strategies. This also holds true for our study. The results suggest that additional a ention to maritime safety management should not be restricted to areas with the highest accident probabilities, but other spatial aspects should also be taken into account. Regional oil spill contingency plans should also take into account the spatial distribution of risk, which might mean, for example, allocating more oil spill-combating equipment to areas with large numbers of highrisk habitats and species.
Furthermore, our fi ndings suggest that this additional risk should be taken into account when planning management actions for specifi c habitats or species. If we want to keep the probability of losing certain habitat types or species small enough, we might need to apply more strict management actions than in situations where there is no additional risk posed by oil spills. This holds true especially if we are not able to control the risk. Some species and habitats (e.g., sandy beaches) are mainly found in areas which are challenging when it comes to combating oil spills, that is, they are located in the outer archipelago or otherwise exposed areas, which may be diffi cult to protect with oil booms especially in bad weather conditions. In these cases it is important to ensure that there are also viable and well-protected occurrences elsewhere in the area, and possibly strengthen their protection with more rigorous regulations.
Our work adds to a broad range of oil spill risk assessment methodologies in the literature, related to both the probability of oiling and its consequences, and takes a step forward toward a comprehensive spatial oil spill risk assessment. With regard to the probability of oiling, the strength of our approach lies in the large number of separate oil spill simulations as well as the probabilistic integration of the scenarios. Although some of the sophisticated oil spill models can be run, for example, with varying weather parameters (e.g., SIMAP;French McCay 2003, modeling is usually still based on predefi ned scenarios thus excluding the possibility of high uncertainties related to the topic. There are some studies that, similar to ours, produce maps describing the probability of oiling based either on historical weather data (e.g., Guillen et al. 2004, Olita et al. 2012, Romero et al. 2013 or a selection of some typical weather scenarios (den Boer et al. 2014 ). However, usually only a few spill locations, spill sizes, or oil types are considered (e.g., Helle et al. 2011, Olita et al. 2012, Romero et al. 2013, Singkran 2013. Further, even if, for example, multiple spill locations are included in the analyses, it is typical to assume that the accident probability is constant across the sites (but see den Boer et al. 2014 ). Our work is a step forward from this, as we have included a BN that off ers an explicit way to handle uncertainty related to oil accidents by describing the variables using probability distributions. Thus, our analysis is not restricted to specifi c scenarios, nor do we have to make any assumptions that all the alternatives are equally probable either.
Regarding the consequences of spills, we have concentrated on threatened species and habitats, the values of which in spatial assessment are based on conservation and legal statuses. Most approaches found in the literature use the www.esajournals.org

HELLE ET AL.
Environment Sensitivity Index (ESI) or a similar approach. In ESI, the coastline is mapped based on the physical properties of the coastline, biological resources sensitive to oil, and the human activities that may be impacted by oil spills, but only the physical elements are ranked quantitatively (Jensen et al. 1998, Petersen et al. 2002. Some studies have given quantitative values and weights also to biological and socioeconomic elements thus producing comprehensive vulnerability indexes (Wirtz and Liu 2006, Castanedo et al. 2009, Fa al et al. 2010, Olita et al. 2012, Frazão Santos et al. 2013. We decided to concentrate on threatened species and habitat types, as they will probably suff er from delayed and uncertain recovery compared to common species and habitats (Ihaksi et al. 2011 ).
Although there are a number of methods that can be used to describe both sides of oil spill risk, that is, the probability and consequences of oiling, approaches that combine the two aspects are still scarce. Olita et al. ( 2012 ) combined a 2yr spill simulation data set describing the probability of oil of hi ing the coast with sensitivity maps, which were based on the geomorphology and the level of coastline protection, and Romero et al. ( 2013 ) used a 30-d oil spill simulation and ESI ranking to produce an index that combines both aspects. Our approach share similarities with these studies, but widens the approach further by also taking into account uncertainties related to accidents.
Yet, there are also some ma ers that need to be taken into account. First, we used a relatively small number of potential accident locations, and increasing them might give us a more comprehensive picture of oil spill risk. However, the current locations are already fairly evenly distributed within the GoF and the AS, and they include relevant locations given our knowledge of maritime traffi c densities as well as species and habitat occurrences. Second, winter is excluded from the analysis, as oil spill models that can reliably mimic the movement of spilled oil in varying ice conditions within the study area are lacking. Winter is, however, a less active period in nature, and thus spills happening in this season might not have as severe negative impacts on coastal ecosystems as spills in other seasons, especially if shoreline clean-up is carried out before, for example, migratory birds start nesting or the growing period begins. If oil spill probability maps for winter become available in future, their inclusion in the modeling framework will be straightforward.
Third, the data used in the study have some limitations. The grid used in oil spill simulations is fairly coarse, and the probability of oiling indicates oil somewhere within a cell, not necessarily in the exact location of threatened species or habitat patches. However, for example, Wikelski et al. ( 2002 ) have shown that even a very low amount of residual oil contamination can increase mortality in wildlife populations a er an oil spill, and we decided to take a conservative approach and assume that oil contamination in the close vicinity of species and habitats is also harmful. Furthermore, we recognize that the database of threatened species is not complete, and there may be occurrences that are not yet, or may never be, known. As the length of the shoreline in the study area is 26 921 km (computed from the national basic map of Finland (1:20 000)), it is unrealistic to assume that the entire area could be covered with extensive surveys in the near future. We argue that somewhat incomplete data sets are not a good reason to refrain from conducting risk assessments, but, as oil spill risk assessment together with risk management are continuous and adaptive processes (Frazão Santos et al. 2013 b ), assessments need to be supplemented with updated data when possible. Also species distribution modeling (Elith and Leathwick 2009 ) might off er an option to improve the ecological data used in the analysis.
To conclude, we have applied a novel, probabilistic method to conduct an oil spill risk assessment for threatened species and habitat types in the northern Baltic Sea. The results show that taking the spatial distribution of biodiversity into account aff ects, for example, our understanding of the most risk-prone areas, that is, a high accident frequency does not necessarily mean higher biodiversity losses and vice versa. Our work, thus, highlights the importance of a thorough risk assessment, which is not based solely on one or two specifi c factors such as accident probabilities or the transport trajectories of spilled oil, but contains also the consequences of spills. We also emphasize the need to take the threat posed by oil spills into account when planning management actions for certain habitat types and species, and our results suggest that www.esajournals.org HELLE ET AL. the risk is unevenly distributed between them. We believe that these fi ndings will be of interest to people around the world working in the challenging fi eld of oil spill risk assessment and management.