Lagrangian Modeling of Arctic Ocean Circulation Pathways: Impact of Advection on Spread of Pollutants Journal of Geophysical Research: Oceans

Sea-ice-free summers are projected to become a prominent feature of the Arctic environment in the coming decades. From a shipping perspective, this means larger areas of open water in the summer, thinner and less compact ice all year round, and longer operating seasons. Therefore, the possibility for easier navigation along trans-Arctic shipping routes arises. The Northern Sea Route (NSR) is one trans-Arctic route, and it offers a potential 10 day shortcut between Western Europe and the Far East. More ships transiting the NSR means an increased risk of an accident, and associated oil spill, occurring. Previous research suggests that current infrastructure is insufﬁcient for increased shipping. Therefore, should an oil spill occur, the window for a successful clean-up will be short. In the event of a failed recovery, the long-term fate of the unrecovered pollutants must be considered, at least until the next melt season when it could become accessible again. Here we investigate the role of oceanic advection in determining the long-term fate of Arctic pollutants using a high-resolution ocean model along with Lagrangian particle-tracking to simulate the spread of pollutants. The resulting ‘‘advective footprints’’ of pollutants are proposed as an informative metric for analyzing such experiments. We characterize the circulation along different parts of the NSR, deﬁning three main regions in the Eurasian Arctic, and relate the distinctive circulation pathways of each to the long-term fate of spilled oil. We conclude that a detailed understanding of ocean circulation is critical for determining the long-term fate of Arctic pollutants. spilled into the ocean. Because of the harsh Arctic environment, an oil spill may not be successfully recovered, so we need to consider where it will go in the following months and years. We released virtual ‘‘particles’’ into a com-puter model of the ocean and tracked their progress for 2 years. In this time, particles traveled, on average, 1,223 km. This demonstrates that pan-Arctic modeling is needed in the event of an unrecovered pollutant spill. Unrecovered oil from one season may be accessible the next spring. By analyzing the spread of our particles, we found that on average 676,917 km 2 would need to be searched to ﬁnd it, but that this is highly dependent on where the spill occurs. Finally, we noted that in some places, particularly the Barents Sea, there was a risk that spilled pollutants could become entrained into deep water, rendering them irrecoverable. particles (a) 9 time a and the and (b) after 2 to assess the fate of irrecoverable oil. Gyre are present in the model results and consistent with the satellite data. This suggests that NEMO accurately simulates the large-scale geostrophic ﬂow at the surface.


Introduction
Marine oil spills are a major concern both environmentally and economically. The financial cost of accidental oil spills can run to billions of dollars, and they have the potential to cause significant damage to marine habitats by contaminating the food web and polluting large stretches of coastline (Carson et al., 2003). It is not possible to completely eliminate the risk of an oil spill occurring, and even thorough clean-up operations can leave some unrecovered oil with environmental impacts (Peterson et al., 2003). Therefore, it is important to understand where spilled oil is likely to be transported to, in order to predict the likely environmental, economic, and social consequences of an unrecovered spill.
An oil spill in the Arctic Ocean is becoming increasingly likely. Permanent Arctic sea-ice is retreating rapidly, and it is predicted that the ocean will be seasonally ice-free by the middle of this century if greenhouse emissions continue at their current rate (Boe et al., 2009;Overland & Wang, 2013;Wang & Overland, 2012). This ''opening up'' of the ocean is fuelling increased interest in using the Arctic for commercial shipping (Aksenov et al., 2017;Lee & Kim, 2015). Figure 1a illustrates one such shipping route, the Northern Sea Route (hereafter NSR), which connects the Atlantic gateway to the Arctic with the Bering Strait and the Pacific Ocean. A schematic of this route is presented in Figure 1.
In turn, this increased shipping activity (Østreng et al., 2013) increases the probability of an accidental oil spill from a commercial tanker or from offshore operations occurring. Winter sea-ice will persist in the Arctic for the foreseeable future, however it is likely to be reduced in thickness and extent, and it will be more mobile (Aksenov et al., 2017). Nonetheless, sea-ice remains a considerable risk for potential shipping accidents, and thus for potential oil spills. The harsh Arctic environment and remoteness of the ocean make this a particularly risky place for a spill to occur.
Once released into the ocean, many factors can govern the fate of spilled oil. These depend on the type of oil released, whether dispersants have been applied, and the environment in which the spill occurs (Afenyo et al., 2015;Mariano et al., 2011). Mixing or dissolution into the water column, photo-oxidation, emulsification, evaporation, sedimentation, biodegradation, and ingestion into the food web are all potential fates for spilled oil (Mariano et al., 2011;Xie et al., 2007).
Oil spills in ice-covered areas behave differently to those in other parts of the world due to their interaction with sea-ice (Afenyo et al., 2015). Ordinarily, oil within the water column, or at the ocean's surface, is transported by ocean currents. However, oil can become trapped and corralled in leads between areas of seaice, and may even become fully ''encapsulated'' into growing ice, effectively isolating it from the ocean below. In this manner, oil can be transported by sea-ice and later, upon ice melt, released back to the ocean in a different location, far from the original spill (Afenyo et al., 2015;Izumiyama et al., 2004).
Leaving aside sea-ice interactions, advection due to ocean currents is the dominant process in determining what will happen to oil that has been mixed, dissolved, or otherwise entrained into the water column in the following months and years (Afenyo et al., 2015). Therefore, understanding the advective pathways in the ocean is key for understanding the long-term fate of spilled oil. The other option is that oil could for at least some of its trajectory be transported by sea-ice, however that is beyond the scope of this work. There is no fixed definition of the NSR, so we have defined a ''main'' route via straits (solid line) and a more poleward ''alternative'' route (dashed). (b) Schematic of Arctic surface circulation. Red: Atlantic inflow following the Arctic Circumpolar Boundary Current (including Barents Sea Branch) and branching at the Lomonosov Ridge. Green: Pacific inflow, following three main pathways: the Alaskan Shelf-break Jet, into the Chukchi Sea through the Herald Canyon, and central flow across the shelf. Note that the flow north of the Canadian Archipelago is currently not definitively established, and presumed flow is shown. (Aksenov et al., 2011). Magenta: Beaufort Gyre. Orange: Transpolar Drift Stream, from Siberia to the Fram Strait.
Biodegradation is a likely eventual fate of an oil spill, but biological processes are inhibited by the freezing temperatures in the Arctic. This means that biodegradation would take longer in the Arctic compared to the rest of the global ocean, meaning that spilled oil would remain an active pollutant for months to years and so long-term consideration of its fate is required (Fingas & Hollebone, 2003). Additionally, due to the short operational season, large distances to ports and other infrastructure, and the generally challenging Arctic environment, there is a significant chance that if a spill occurs, it will not be fully recovered before the winter freezeup makes it inaccessible. During this time, ocean currents and sea-ice can transport oil hundreds of kilometers away from the initial spill location (Main et al., 2016).
Aside from in the Arctic, recent major oil spills have included the Deepwater Horizons incident in the Gulf of Mexico. Ocean models were used to investigate the spread of pollutants here (Macfadyen et al., 2013;North et al., 2013;Weisberg et al., 2017), including Lagrangian analysis of modeled currents. It was found that ocean currents played a dominant role in determining where oil would go, and this held until oil reached the coast, when Stokes drift became the important factor for determining how the oil would beach (Weisberg et al., 2017).
To date, no major spill has occurred in the Arctic. An oil rig ran aground in the Chukchi Sea in 2012, but this did not lead to actual pollution (Meier et al., 2014). However, there has been one particularly notable oil spill in ice-affected waters: In 1989, a major spill from the Exxon Valdez tanker occurred off the coast of Alaska (Peterson et al., 2003). One contributing factor to this event was the tanker deviating from its normal shipping lane to avoid icebergs. This caused an accident, and 11 million gallons (50,000 m 3 ) of Prudhoe Bay crude oil was spilled into the ocean. The financial cost of the spill ran to billions of dollars, with Exxon spending over $2 billion on oil spill response and restoration (Carson et al., 2003).
Consequently, preparedness modeling around ocean circulation in the Arctic is vital to understand where oil spilled in the region will be located in the springtime when it becomes accessible again, both to permit estimation of the level of the potential recovery costs, and to understand the likely domain and severity of environmental and ecological impacts should a spill similar to Exxon Valdez occur within the Arctic Ocean.
By way of a summary, large-scale circulation pathways in the deep Arctic Ocean are primarily driven by the wind and by the inflows from the North Atlantic and Pacific oceans (Aksenov et al., 2011;Pnyushkov et al., 2013Pnyushkov et al., , 2015Proshutinsky et al., 2015;Proshutinsky & Johnson, 1997). Since the ocean stratification below the Arctic halocline is week, the intermediate depth currents are strongly influenced by the by the oceanic ridges and steep topography of the continental shelf slope (Aksenov et al., 2011). The exchanges between the shelf and the deep part of the ocean occurs through cascading of the dense shelf waters Ivanov & Golovin, 2007) and the cross-slope currents, driven by the along-shelf component of the wind stress through the Ekman transport mechanism .
The presence of sea-ice serves to decrease wind forcing and makes the ocean circulation relatively slow, however this ''inhibition'' is anticipated to decline as the ice retreats. In bathymetric terms, the Arctic Ocean is roughly evenly divided between shelf seas (up to 200 m depth) and deep ocean, with the latter split by the Lomonosov Ridge into the Amerasian and Eurasian basins (Bj€ ork et al., 2007).
As Figure 1b illustrates, the anticyclonic Beaufort Gyre (magenta) dominates the circulation in the Amerasian Basin (Bluhm et al., 2015;Proshutinsky & Johnson, 1997). Meanwhile, cyclonic regimes govern circulation in the Eurasian Basin (red), with currents guided by local bathymetry, and following shelf-breaks and ridges (Carmack & Wassmann, 2006;Wassmann et al., 2015). Between the cyclonic and anticyclonic regimes of the two basins, the Trans-Polar Drift Stream (Orange) carries water from Siberia to the Fram Strait across the deep Arctic (Bluhm et al., 2015).
The NSR shown in Figure 1a predominantly crosses the Eurasian Arctic shelf. This region is comprised of five seas: the Barents, Kara, Laptev, East-Siberian, and Chukchi. The Barents and Chukchi seas are inflow shelves, and are influenced by incoming Atlantic and Pacific water, respectively (Carmack & Wassmann, 2006;Williams & Carmack, 2015). Atlantic water entering the Barents Sea Opening flows east across the shallow shelf. Dense water is formed in the Barents Sea, sinks, and becomes a branch of the Arctic Circumpolar Boundary Current (ACBC) (Aksenov et al., 2011). The ACBC has three cores: the main of which is a lower shelf slope current centered at around the temperature maximum in the Atlantic Water (AW) layer at 300 m depth, originated from the AW inflow through Fram Strait, a deeper Barents Sea Branch at ca. 1,000 m, and as well as a surface/subsurface branch (down to 150 m depth in the water column), aka Arctic Shelf-Break Branch (ASBB) (Aksenov et al., 2011). We refer to all these as part of the ACBC. We assert that the above boundary current structure has been observed by the NABOS and other observational programs and has been simulated in high-resolution models.
At the other end of the Eurasian shelf, a smaller amount (1 Sv) of Pacific Water enters the Chukchi Sea via the Bering Strait (green in Figure 1b), where it can follow the Alaskan Shelf-Break Current toward Canada, flow into the East-Siberian Sea and potentially join the transpolar drift downstream, or enter the Beaufort Gyre (Aksenov et al., 2011;Carmack & Wassmann, 2006).
Between the Barents and Chukchi seas lie the interior seas of the Arctic (Kara, Laptev, and East-Siberian). These are influenced by freshwater runoff from the Ob and Yenisei near the coast, and from the ACBC along the shelf edge in the north (Carmack & Wassmann, 2006;Williams & Carmack, 2015). Wind forcing plays an important role here, with cyclonic summer atmospheric circulation favouring an eastward transport in the Laptev Sea. (Bauch et al., 2009). Significantly for communication with the deeper interior of the Arctic Ocean, Ekman pumping drives upwelling and downwelling at the shelf-break and associated cross-shelf exchange (Williams & Carmack, 2015).
While this sketch describes the general pattern of ocean circulation in the Arctic Ocean, the region is unsurprisingly also characterized by strong seasonal and interannual variability (in addition to the secular trend driven by anthropogenic climate change). All of these factors combine to make the Arctic Ocean a complex and variable region for understanding and planning spill responses.
In the following section, we outline the Lagrangian modeling technique used to investigate the impact of advection on a potential Arctic oil spill from various locations along the NSR (marked in Figure 1a). We then quantify this with various metrics, and discus the consequences of ocean circulation for the spread of pollutants from in the Arctic Ocean.

NEMO (Nucleus for European Modeling of the Ocean)
In this study, we use the ORCA0083 1/128 resolution configuration of the NEMO (Nucleus for European Modelling of the Ocean) general circulation model (GCM) coupled to the Louvian-la-Neuve Ice Model (LIM2) sea-ice model (Fichefet & Maqueda, 1997;Goosse & Fichefet, 1999;Madec, 2014). Here we present model description relevant for this study, for more detail the reader is referred to Madec (2014).
NEMO is a global z-level model with a fully nonlinear free surface. Horizontal resolution in the Arctic is 3-5 km, making it eddy-resolving in the central Arctic Ocean but only eddy-permitting on the shelves due to the small Rossby radius of deformation (Nurser & Bacon, 2014). The model has 75 vertical levels, with spacing varying from 1 m at the surface to 204 m at 6,000 m (there are 31 model levels in the upper 200 m). The model simulates vertical mixing using the turbulent kinetic energy (TKE) mixing scheme (Blanke & Delecluse, 1993) and uses the total variance dissipation (TVD) advection scheme for active tracers (Madec, 2014). The LIM2 ice model uses Elastic-Viscous-Plastic rheology (EVP) (Hunke & Dukowicz, 1997), implemented on a C-grid (Bouillon et al., 2009), with thermodynamics based on two layers of ice and one layer of snow (Fichefet & Maqueda, 1997). It is coupled to the ocean model at every ocean time step, with a nonlinear quadratic drag of the sheer between the ice and ocean. The model has been found to accurately simulate sea-ice in the Arctic Ocean (Aksenov et al., 2017;Johnson et al., 2012;Wang et al., 2016).
The model is forced at the surface boundary using the DRAKKAR forcing set (DFS) atmospheric reanalysis (Brodeau et al., 2010). This is comprised of 6 hourly data for atmospheric winds (from the ERA40 reanalysis), temperature and humidity, daily radiative fluxes (shortwave and longwave) and monthly means for precipitation (rain and snow; from the CORE2 reanalysis) and runoff (Brodeau et al., 2010;Timmermann et al., 2005). In the simulation used here, NEMO was run from rest, with forcing from the beginning of 1978 until the end of 2015, and output saved as 5 day means. NEMO is widely used by the research community for global studies at a variety of resolutions, including ORCA0083, (e.g., Duchez et al., 2014;Janout et al., 2015;Marzocchi et al., 2015;Srokosz et al., 2015). In the Arctic, it has been used extensively in the coarser 1/48 configuration (e.g., Aksenov et al., 2017;Lique et al., 2010;Popova et al., 2013). Evaluation of the circulation in NEMO (ORCA025), by way of calculating the barotropic stream function across the Arctic is presented in Lique et al. (2010), though the authors note the difficulty in accurately observing surface currents as sea surface height cannot be directly observed (Lique et al., 2010)-although reanalysis products are available and used here.
Although the recent realization of ORCA0083-N06 used here has undergone extensive validation globally, it has not been comprehensively evaluated in the Arctic. Further evaluation of the 1/128 NEMO ORCA0083-N06 run is presented in section 3.1. Here we compare ice cover against satellite-derived reanalysis data. We also compare modeled sea surface height with satellite observations. This is used in conjunction with an analysis of the model's barotropic stream function to evaluate the simulated circulation.
As with all models, NEMO is not without its limitations. For example, the configuration of NEMO used in these experiments lacks tides and wave model. We do not expect these to have a significant effect on our results, but it does place restrictions on the conclusions that we can draw in some coastal areas where tides play a particularly significant role in the ocean dynamics (Luneva et al., 2015;Padman & Erofeeva, 2004). Specifically, tides and waves are both important for mixing oil into the water column and dispersing it. However, given that we are modeling advective pathways-rather than directly modeling the physics of the oil itself, which would vary strongly based on the type of oil spilled-mixing is not directly accounted for anyway. We instead aim to describe the circulation pathways followed by the ocean currents in order to give a more general overview of where spilled oil could go, assuming it has already become mixed into the water column.

Lagrangian Modeling
There are two approaches by which ocean models can address pollutant dispersal: (1) online representation via a passive tracer, whose concentration is determined by the resolved circulation and the parameterized mixing-for more detail, see Madec (2014); and (2) offline, using saved output from a preexisting run of the model. This approach uses Lagrangian ''particles,'' whose positions in space and time are updated according to the saved mean velocities and does not require the full model to be rerun (Blanke & Raynaud, 1997). Both have advantages and drawbacks.
The transport of online tracers is consistent with model transport processes of advection and diffusion, and mixed layer processes such as convection. NEMO employs Eulerian meshes to solve the tracer evolution equation numerically, thus this approach is often called ''Eulerian.'' However, this Eulerian approach comes at a significant computational cost because it requires rerunning the high-resolution global model itself for each simulated spill scenario.
In contrast, by making use of output from an already existing (and computationally expensive) simulation, and by calculating for only a fixed number of trajectories, offline Lagrangian particles require significantly less computational resource. This reduced cost can particularly suit studies where trajectories are repetitively initialized from multiple sites at multiple time points, a situation which would require many separate simulations for online tracers. Using hundreds of thousands of discrete trajectories, it is possible using Lagrangian particles to identify advective pathways and their variations. The downsides of using this offline approach are that individual Lagrangian particles effectively represent large quantities of pollutant, and parameterized mixing processes are unrepresented, so we can only consider the effects of advection but not diffusion. This can be compensated for by using many particles over an ensemble of releases. Furthermore, while the Lagrangian method does allow for subduction due to nonzero vertical velocities, it cannot account for convection-i.e., vertical mixing due to water-column instability. This is most likely to affect the results in the areas of Atlantic inflow where convection is prevalent.
For the purposes of this study, we use the ARIANE particle-tracking software package which uses the Lagrangian methods outlined above to calculate the evolution of trajectories (Blanke & Raynaud, 1997). Here ARIANE reads in 3-D velocity fields from the 5 day mean NEMO output, and uses this to disperse virtual ''particles'' released into the model's flow field. These particles are transported per a bilinear interpolation of the velocity field, using an analytical method to solve for particle translation through model grid cells, and the resulting trajectories are stored for analysis at daily frequency. Although ARIANE does not include horizontal mixing (which, in part, represents subgrid processes in the ocean model), the high resolution of NEMO used here accounts for most of the relevant transport explicitly. ARIANE is mass-conserving and powerful for describing the large-scale, long-term impact of advection. However, it does not account for turbulent mixing, and since it works with 5 day mean advection fields, it does not guarantee that particles will exactly follow constant density surfaces.
Individual trajectories from each simulated pollution event can be plotted to highlight the different pathways followed by each particle, with the distribution of particles after some given time indicating dispersal size, distance, and shape. Using multiple releases initially close in space and time can provide information on the uncertainty associated with different spill sites and dispersal routes. Here we use the term ''advective footprints'' to describe the ensemble of trajectories from a given release site. Lagrangian analysis, specifically using ARIANE and the advective footprint approach advocated here, has previously been employed to model oil spills and other pollutants (Main et al., 2016;Robinson et al., 2017).

Experiment Design
In order to evaluate the impact of advection of a potential Arctic oil spill, we consider the parts of the ocean most at risk of an accident occurring. The Northern Sea Route (NSR) is a major shipping corridor running between the Barents Sea in the west and the Bering Strait in the east (Lee & Kim, 2015). It offers a connection between North West Europe and the Far East which, in the future, could be more economically viable than the normal shipping route via the Suez Canal (Liu & Kronbak, 2010;Schøyen & Bråthen, 2011). It offers a 40% reduction in distance compared to the Suez route (Schøyen & Bråthen, 2011), and this has the potential to cut sailing times by up to 10 days (Aksenov et al., 2017). Both the reduced time at sea, and associated saving in fuel, make this a potentially attractive route for shipping companies.
To model the impact of advection on an oil spill from the NSR, virtual Lagrangian particles were released into the NEMO flow field at 15 sites along the Northern Sea Route (see Figure 1a). These locations were selected to sample the length of the NSR, between Murmansk in the west and the Bering Strait in the east, via the Barents, Kara, Laptev, East-Siberian, and Chukchi Seas. The NSR is not strictly defined as a single corridor, and different definitions have been used in the previous literature (e.g., Aksenov et al., 2017;Lee & Kim, 2015). Here we take this into account by defining a ''main'' route via straits connecting each sea, as well as a more northern ''alternate'' route.
Release sites at Murmansk and the Bering Strait were selected as the start and end points of both routes. The straits along the ''main'' route are potentially high-risk locations due to variable sea-ice cover, shallow water and restricted room for manoeuvre, which could increase the risk of a shipping accident. A further five sites were taken at the midpoints between these straits to sample each of the Eurasian Arctic seas. Finally, the borders between each sea were sampled from our northward ''alternate'' route to complete the set of release locations. The locations of each release site are marked in Figure 1a.
At each of the 15 sites, 100 particles were seeded over a regular 10 km 3 10 km grid, chosen to represent an area covered by oil after several days of spreading. These releases were repeated every 10 days throughout the navigable season of the NSR (taken to be start of June to end of October), every year from 1990 to 2009, to take into account seasonal and interannual variability of the circulation. All releases from a given year were conducted as part of the same experiment, resulting in twenty experiments of 22,500 particles.
All trajectories were advected with the full 3-D velocity field for 2 years from their respective launches, with the particles' positions recorded daily. We consider this 2 year period to account for spills which are not recovered before the Arctic freezeup begins, and therefore persist in the Arctic for beyond one winter. Previous research has evaluated dispersal over shorter time periods (1 year) but this has been in regions with warmer waters and no sea-ice (Main et al., 2016).
It is important to note that these particles do not directly represent real particles of oil. With only 100 particles per release, each is potentially representative of a large quantity of pollutant (see section 3.3.6. on the sensitivity of this assumption). Also, the particles are modeled here as being neutrally buoyant, and represent oil that is dissolved (or emulsified) within seawater rather than that which is floating as a sea surface slick or whose density changes with time.

Transport Metrics and Advective Footprints
We begin by simply plotting trajectories from each release site in order to give a qualitative description of the pathways that they follow (see section 3.2.1). Additionally, we define quantitative metrics to describe the distances traveled by particles, the uncertainty associated with where they go, and the sensitivity to when a spill occurs: First, we are interested in how far particles go. This is trivial to quantify, and we used two metrics: total distance traveled (the sum of all the distances traveled in each time step) and A-B distance traveled by particles from each release site. These distances were calculated after (a) 9 months, to represent the typical time between a spill occurring in the summer and the melt season beginning next spring, and (b) after 2 years, to assess the fate of irrecoverable oil.
In order to better describe the spread of Lagrangian particles, we introduce the concept of ''advective footprints.'' Footprints are defined as the area of ocean covered by Lagrangian particles some time T after release. We choose T 5 9 months (270 days) to correspond the approximate time between an oil spill occurring and the beginning of the melt season next spring.
To calculate the area covered by particles, we define a grid and count the number of cells occupied during the 270th time step (i.e., only cells occupied in the 270th time step, and not counting cells passed through on their way). We use a coarsened version of the ORCA12 grid-10 cells in the i direction and 11 in the j direction. This was chosen to provide approximately regularly sized square cells, and of a size that typically produced continuously filled footprints. The ''area of advective footprint'' is then defined to be the sum of the areas of each occupied grid cell. This is calculated for each release (so 15 releases for each of the 20 years studied) individually. This figure is then averaged over the 300 releases for each site to give an estimate of the area that would likely be affected by a spill from a particular location-i.e., a measure of the horizontal spreading of our particles. We then compare how this value varies with respect to both season and year of launch to assess the interannual and intraannual variability of our experiments.
Additionally, we consider the ''envelope'' of all of these footprints-i.e., the total area covered by all 30,000 particles from a given launch site. The size of this ''envelope'' is representative of the uncertainty associated with pathways from a given location. Spilled oil could go anywhere reached by the Lagrangian particles, but it will not necessarily follow every pathway: each trajectory represents one of many possible pathways. This metric gives an overview of the areas potentially at risk.
Finally, we investigate the likelihood of subduction and the pathways followed by subducted particles. Oil that is entrained deep into the water column would probably not be recoverable, therefore it is necessary to understand where there is an enhanced risk of this happening. We take 100 m as an approximation for the maximum mixed layer depth in the Arctic Ocean, and compare the fraction of particles below this threshold after 9 months. While this is deeper than the shallow mixed layer depth in the Arctic Ocean, this threshold was chosen to ensure that particles reaching these depths were definitely below the upper mixed layer regardless of where in the Arctic they sank below 100 m.

Model Evaluation
The global NEMO model has previously been validated in 18, 1/48, and 1/128 resolution configurations (Marzocchi et al., 2015). Here we use the latest (as of 2017) run of 1/128 NEMO, although thorough evaluation of the Arctic at this resolution has only been done in the Kara and Laptev seas . However, Arctic circulation and exchanges have been extensively validated in the coarser 1/48 version, with the model reproducing observed features Lique et al., 2010;Popova et al., 2013). Nonetheless, model evaluation is an ongoing endeavor, and here we focus on the performance of key model features that are most pertinent to the present study: ocean circulation and sea-ice cover.

Ocean Circulation
In order to evaluate the surface currents produced by the NEMO model, two variables are investigated: sea surface height (SSH) and the barotropic stream function (BSF).
Sea surface height is verified against satellite Dynamical Ocean Topography (DOT) inferred by the CPOM (Centre for Polar Observation and Modelling) Envisat (2003-2011) (Armitage et al., 2016). The SSH field is used as a pan-Arctic proxy for observed and modeled surface geostrophic currents. To compensate for the two data sets being measured relative to different reference level-a geoid for satellite data and a spherical reference SSH for the model-the anomalies rather than absolute SSH are compared and presented in Figure 2.
In both the satellite-derived and modeled data, the magnitude of the differences in sea surface height (SSH) between the highest and lowest parts of the Arctic is approximately 1 m on average between 2003 and 2011. Key features such as the raised Beaufort Gyre are present in the model results and consistent with the satellite data. This suggests that NEMO accurately simulates the large-scale geostrophic flow at the surface.
The barotropic stream function was calculated for the average NEMO velocity fields (not shown, though the barotropic stream function for the 1/48 version is discussed thoroughly in other work, Lique et al., 2010). It shows a similar spatial pattern to the SSH for the same period. Because most of the stratification in the Arctic is confined to the upper 200 m, the ocean circulation at intermediate depth is largely barotropic (Aksenov et al., 2011;Pnyushkov et al., 2015), and supports the interpretation that the barotropic flow in the model is also qualitatively in agreement with the satellite-inferred data-see Figure 2.

Sea-Ice
Further evaluation is provided by comparing modeled and reanalysis (from the National Snow and Ice Data Center (NSIDC)) sea-ice extent. This is presented in Figure 3. Two measures are used to analyze the horizontal extent of sea-ice cover: we compare the mean summer (June-September) ice extent in the NEMO and  NSIDC data sets (Figure 3a), and we look at how well NEMO reproduces the seasonal cycle of ice growth/ decline (Figure 3b). Both metrics are compared between 1900 and 2012 to cover the period of our Lagrangian experiments. Sea-ice extent is taken to be the area covered with ice concentration greater than 15%.
Generally, NEMO reproduces real ice conditions well. The current downward trend in summer sea-ice extent is clearly visible (Figure 3a) with minima in 2007 and 2012 minima reproduced. Figure 3b shows the monthly cycle of ice growth and decline, were NEMO is consistently within one standard deviation of the reanalysis data.
NEMO's seasonal cycle is also in good agreement with the NSIDC reanalysis data (Figure 3b). The total extent of the modeled ice cover is typically an overestimate, especially during the winter months, but it is accurate to within two standard deviations of NSIDC reanalysis data. Spatially, the ice cover in NEMO and NSIDC is reasonably similar, with sea-ice present/absent in approximately the same areas.
For a more detailed evaluation of NEMO/LIM2, the reader is referred to Wang et al. (2016).

Lagrangian Experiments 3.2.1. Arctic Circulation Pathways and Time Scales
In order to assess the advective pathways likely to be important for an Arctic oil spill, Lagrangian particles were released from 15 sites along the Northern Sea Route. Releases were performed every 10 days from May to October, from 1990 to 2009. Examples of trajectories from each site are presented in Figure 4, showing all releases from 2000. Note that reds (corresponding to the earlier parts of trajectories) are plotted over yellows/greens/blues from later in the particles' journeys to highlight the fastest advective time scales.
From Figure 4, the first result that is immediately apparent is that, depending on where the particles were released 1 year can be sufficient to reach the center of the Arctic Ocean. Within 2 years it is possible for them to leave the Arctic completely and flow into the North Atlantic. The spatial scale of these advective pathways confirms that a pan-Arctic outlook is necessary to understand the long-term fate of potential oil spills from the NSR.
As is also clear from Figure 4, advective pathways in the Arctic are highly variable depending on where (and when) Lagrangian ''particles'' are released. First, we address the sensitivity to release site. Below we detail the spread of the particles through different regions of the Arctic Ocean and Arctic seas: 3.2.1.

Barents and Kara Releases
The north-eastward flow of the Atlantic Water through the Barents Sea, the Barents Sea Branch (Aksenov et al., 2010(Aksenov et al., , 2011, is the main circulation feature, and Lagrangian particles released in the Barents Sea show this pathway (see Figures 4a and 4b). The particles followed two main trajectories; (1) eastward flow within the Murmansk Current, then through the Kara Gate into the Kara Sea, and (2) flow northward of Novaya Zemlya toward the Eurasian shelf-break within the West Novaya Zemlya Current, where they flow via the St. Anna Trough and across the northern Barents Sea shelf into the deep Arctic Ocean (Aksenov et al., 2010). Trajectories are relatively fast in this part of the Arctic, with 4-5 month transient times to cross from Murmansk to the Kara Gate. The trajectories also map the cyclonic gyre in the Central Basin of the Barents Sea, with some of the particles entrained in the gyre and following the westward flow toward the Barents Sea Opening (cf. Figure 4a and Aksenov et al., 2010).
As in the Barents, northward and eastward currents dominate in the Kara Sea (Figures 4c-4e). The particles released in the Kara Sea travel via the St. Anna Trough in the Nansen and Amundsen basins of the Arctic Ocean (Dmitrenko et al., 2015), and also via the Vilkitsky Strait Current into the Vilkitsky Strait and the western Laptev Sea . There is an episodic flow back from the Kara Sea through the Kara Gate into the Barents Sea. This is the Litke Current, observed on some occasions (Pfirman et al., 1997), and found to be correlated with an atmospheric sea level pressure gradient across Novaya Zemlya between the Barents and Kara Seas, and by seasonal variations in the buoyant eastward current at the southern end of the Kara Gate.
After leaving the St. Anna Trough, particles go on to follow the Eurasian shelf-break within the Shelf-Break Current (SBC). This is one of the cores of eastward flow, guided by topography Arctic Circumpolar Boundary Current (ACBC; Aksenov et al., 2011), and along the rest of their trajectories subduction down the shelf slope was seen to occur (Figures 4c-4e). From there, particles recirculate in the Nansen and Amundsen basins (Figures 4e and 4f), eventually heading to exit the Arctic Ocean via the western part of Fram Strait toward the Atlantic. Very few particles released in the Kara Sea reach the Atlantic within the 2 years studied.

Laptev, East-Siberian and Chukchi Releases
Lagrangian particles released on the Laptev Sea shelf depict the well-established circulation pattern: they flow toward the shelf slope, where they can follow one of two routes: the ACBC or the Transpolar Drift Current (hereafter Transpolar Drift) from Siberia to the Fram Strait (Aksenov et al., 2011;Dmitrenko et al., 2005;Janout et al., 2015;Newton et al., 2008) (see Figures 4f-4i). Trajectories following the ACBC are strongly guided by the seabed topography of the Arctic, and branch away from the shelf-break at the Lomonosov Ridge (Woodgate et al., 2001). Here they join the Transpolar Drift and flow toward the Atlantic. Particles released from all sites around the Laptev Sea are able to reach the Atlantic Ocean within 2 years.
Trajectories in the East-Siberian Sea are notably slower than those in other shelf seas (see Figure 4l in particular). Currents in the area vary in time and space, and the particles are carried in all directions, owing to oppositely directed inflowing waters from the Laptev Sea on one side (which carries particles east when it dominates) and the Chukchi on the other (which carries particles west) (Timmermans et al., 2014). Particles from the East-Siberian Sea remain on the shelf for a considerable time, typically beyond a year. Once these trajectories enter the central Arctic Basin, the Transpolar Drift is the principal pathway. Circulation in the Chukchi Sea is driven by inflow from the Pacific and by the wind Timmermans et al., 2014). The majority of particles in our experiments follow Pacific Water pathways, of which three main routes are noted: (1) the western route into the Chukchi Sea through the Herald Canyon, (2) the eastern route via the shelf-break jet through the Barrow Canyon and along the Alaskan Shelf-break (Appen & Pickart, 2012) and (3) the route across the shelf crossing the Herald Shoal and the Hanna Shoal via the Central Channel (Timmermans et al., 2014). Trajectories originating in the Chukchi Sea can become entrained into the Beaufort Gyre within 2 years. This is illustrated in Figures 4m-4o. Further detail of the variability of these pathways is discussed in section 3.2.5.

Distance Traveled
As a first-order metric, we are interested in how far particles go. The mean distances (across all experiments) traveled by particles from each launch site is presented below in Table 1. The distances traveled after 9 months and after 2 years are presented, both as direct A-B and the total length of the path traveled for each of the 15 release sites used, listed from west to east. The full path is a proxy for average current speed along a trajectory, whereas the A-B distance represents how far away from the initial spill site a recovery operation would need to consider.
From Table 1, the distance traveled shows only modest sensitivity release site when the full path length is considered. All sites fall within 300 km of the mean distance traveled (1,497 km) across all sites for the first 9 months of advection: i.e., particles are transported at approximately the same speed regardless where they are seeded from.
However, this is not the case when considering the direct A to B distance (calculated as the shortest distance along the surface of the Fischer Spheroid) for each site. Excluding the Bering Strait, particles starting from release sites at the west of the Northern Sea Route traveled significantly further (typically about twice as far) as those from the eastern end of the NSR (excluding the Bering Strait site, which was more in line with the western sites).
Here the cutoff between west and east was the boundary between the Laptev and East-Siberian Seas. Particles originating from the Barents, Kara and Laptev seas traveled on average 771 km in nine months and 1,454 km in 2 years, whereas those from the East-Siberian Sea (including New Siberian Islands and Sannikov Strait release sites) and Chukchi sea (excluding Baring Strait site) traveled an average of 382 km in nine months and 834 km in 2 years.  The marked difference in A-B distance despite the similar total path distance demonstrates that particles launched toward the west of the route (and also in the Bering Strait) follow more direct routes, whereas those in the east are more prone to recirculation.

Advective Footprints: Horizontal Spread and Uncertainty
We now characterize the horizontal spread of particles. This is done by dividing the ocean into a grid (full details in section 2.5) and counting the total area of occupied grid cells. Here we consider grid cells that are occupied 9 months after particles were initially released. Again, 9 months was chosen to represent a typical time between a spill occurring and the next melt season beginning. In the context of an oil spill, this corresponds to the area that should be expected to be reached by potentially contaminated waters.
Each release of 100 particles per release site is considered separately (so a total of 300 releases per site: 20 years of releases with 15 releases per year) and these are averaged to give a typical release footprint size for each release site. Particles were initially distributed over an area of 100 km 2 . These footprint sizes are presented below in Table 2. Table 2 shows that the spreading of particles is strongly dependent on where they are initially released from. There is a clear east/west divide, with western release sites (1-5) typically producing larger footprints and central (6)(7)(8)(9)(10)(11) and eastern (12)(13)(14)(15) ones producing smaller footprints.
The smallest footprints were associated from the middle of the Northern Sea Route, around Siberia in the Laptev and East-Siberian Seas (including the New Siberian Islands and Sannikov Strait release sites between them). The East-Siberian Sea site produced anomalously small footprints, on average less than 10,000 km 2 .
The average area for the other three sites in this group was 27,267 km 2 (average standard deviation 15,518 km 2 ), making them only half the size of the average footprint area for the remaining 11 release sites.
The five westernmost release sites (Murmansk-Kara Sea) produced the five largest footprints. The mean area of footprints from these sites was 67,866 km 2 (average standard deviation 10,665 km 2 ); 48% bigger than the average for the other release sites, excluding the anomalously small footprint from the East-Siberian Sea release site. Including all sites between Murmansk and the Vilkitsky Strait, the average area of footprints from the west of the NSR was 61,361 km 2 (average standard deviation 12,317 km 2 ).
Considering each of the 300 releases from each of the 15 releases sites, the average area of advective footprints 9 months after particles were releases was 45,762 km 2 (average standard deviation 13,307 km 2 ). Excluding the five largest and four smallest sites highlighted, the typical area of a footprint was 44,161 km 2 (average standard deviation 15,233 km 2 ).  Additionally, the area of the ''envelope'' surrounding all 300 footprints from each site (again after 9 months) was calculated, i.e., the total area affected by at least one release. While the first metric addresses how waters from a particular release spread, this metric addresses part of the uncertainty in where they will go. As with the individual footprints, western release sites were found to have the largest ''envelopes'' while the smallest were in and around the East-Siberian Sea. The large footprint ''envelopes'' in the west were typically associated with particles being rapidly entrained into a well-organized current (the ACBC), which enabled them to travel further (see Table 1) and hence produce a larger envelope than particles which stay closer to their release sites.
Finally, we calculate the ratio of the ''individual footprints'' and the ''envelopes.'' This demonstrates that the uncertainty in where spilled oil could go (the envelope) is an order of magnitude larger than the area over which contaminated waters would be expected to spread (the individual footprints). This ratio provides a measure of the variability of advective pathways from a given site: larger ratios correspond to footprints with less overlap between experiments, suggesting more variable circulation. The combination of the especially large envelope and ratio from the Severnaya Zemlya site is partly due to the highly variable circulation around this area-discussed further in section 3.2.5. The largest ratio was produced by the East-Siberian Sea releases, but it should be noted that this was due to both an abnormally small envelope and individual footprints, which stay relatively close to their initial release site.

Subduction
Having addressed the horizontal spread of our particles, we now look in the vertical direction. It is important to note that we are only considering subduction due to advection, as our Lagrangian technique does not explicitly include convection. All particles were initially released at the ocean surface, but they were not constrained to stay there. Here we investigate the fraction of particles which sink below 100 m depth. This was chosen as a rough approximation for the depth of the upper mixed layer, as any oil entrained deeper into the water column will be especially difficult to recover/unrecoverable. Two snapshots were investigated: the time step after 9 months of advection and the time step after 2 years of advection, as with the distance traveled metrics. At these snapshots, particles were classified as either above or below a 100 m threshold. The results are presented below in Table 3.
From Table 3, it is clear that subduction is not a major concern for the majority of releases locations, with two notable exceptions: the Barents Sea and Novaya Zemlya sites (and, to a lesser extent, Murmansk and the Bering Strait). One in five particles seeded from the first two sites ends up below the 100 m threshold after 9 months of advection, and is potentially unrecoverable.
No particles from any release in any year were below 100 m after 9 months from four central release sites: these were in the Laptev and East-Siberian Seas, as well as the two release points on the border between these seas. We then investigated where this subduction occurs. An illustrative example, with trajectories colored to indicate their depth, is presented in Figure 5. In this figure, trajectories from 2007 are shown in three different groups, to highlight the regional differences in subduction.
From Figure 5, we can see that subduction mainly occurs for the western group of release sites, as noted in Table 3. It is apparent that this occurs mostly around the Eurasian shelf-break, where particles downwell across the shelf slope. These deep trajectories are guided by the bathymetry of the Arctic Ocean and tend to follow the Arctic Circumpolar Boundary current along the Eurasian shelf-break.

Sensitivity to Time of Release
The distance traveled, ''advective footprints'' and subduction metrics discussed in the previous sections were compared with respect to launch site. Here we repeat that analysis, but instead of comparing different release locations, we compare how these metrics vary with respect to time of release. This is presented in Table 4: first we compare the sensitivity to year of release (averaging over all release sites Note. The history of trajectories is not considered here, only whether they are above or below the threshold for the two snapshots studied. and all releases within the given year) and then across each of the 15 releases per year (averaging over all release sites and all years). For each of the three metrics, we look at a snapshot 9 months after particles were initially released.
From Table 4, it is clear that none of these metrics show significant sensitivity to year of release. Interannual variability is present, but with no clear trend in any of the metrics used here.
However, two of the metrics showed significant intraannual variability. The mean footprint area metric suggests slightly reduced spreading of particles correlating with later releases. The most notable trend came in the subduction metric. Here the fraction of particles subducted below 100 m steadily increases toward the end of the season, and it more than double between the first set of releases in June and the last set of releases in October.
As noted in section 3.2.5, subduction is not a major occurrence for most release sites. However, for those where it does matter, it is more likely to happen for particles launched in the autumn than those starting in the spring.
Investigating the sensitivity to time of release highlighted variability in the pathways that trajectories follow. Multiple major pathways exist for some sites, and three examples (specifically the Bering Strait, Kara Gate, and Severnaya Zemlya releases) with clear contrasts are presented in Figure 6.  Figures 6a and 6b show particles which were only launched 20 days apart from each other, yet flow in markedly different directions. Figure 6a highlights the Alaskan Shelf-Break current, which leads to trajectories along the North American coastline, before seeming to join the Beaufort Gyre. Figure 6b shows the other extreme, with particles flowing into the Chukchi Sea and affecting the Russian coastline. Figures 6c and 6d show a more subtle variation. As previously noted, the dominant flow in the Kara Gate is from the Barents into the Kara Sea, but occasionally this is reversed. An example (Figure 6c) which shows this ''reversed'' flow west through the Kara Gate and then north-eastward along the Novaya Zemlya coast is shown. Figure 6d shows the usual eastward-only flow. This could be indicative of a wind-driven blocking event, analogous to similar events in the Vilkitsky Strait . Further investigation (not shown) suggests that the usually dominant buoyancy driven current (west to east) interplays with a wind-driven current in the opposite direction, which correlates with increased atmospheric pressure to the north-west of Novaya Zemlya.  Figure 6e, and all bar one follow the Arctic Circumpolar Boundary Current in Figure 6f. Both pathways are seen frequently, as this marks a region of variable flow: it is the boundary between where the eastward flow dominates in the Barents and Kara seas, and the east where the Transpolar Drift Stream is more common. 3.2.6. Sensitivity to Initial Area of Release So far in this study, we have used 100 km 2 as the initial area covered by our Lagrangian ''particles.'' This was chosen as a first-order approximation for the size of a typical oil spill. We have investigated the spread of these particles and how this varies depending on when and where particles are released, but it is also important to assess how our choice of initial conditions could affect these results. In order to do this, we repeat our experiment for a typical year (2000 was chosen as it is in the middle of the period studied), but with particles released over different sized areas: 400 km 2 (20 km 3 20 km), and 25 km 2 (5 km 3 5 km). We keep the number of trajectories the same in each case. Exactly as before, we then calculate the area covered by particles 9 months after their release. We also compare this to the results from year 2000 in the standard 10 km 3 10 km initial release experiments.
Differences in dispersion were apparent but small, with larger release areas leading to slightly larger distributions after nine months, see Table 5. Results from the 5 km 3 5 km, 10 km 3 10 km, and 20 km 3 20 km experiments are compared: It is apparent from Table 5 that the sensitivity to size of initial release is low. Excluding the East-Siberian Sea release site, which had an anomalously low footprint area anyway, the experiments from the largest release boxes spread to an area only 30% larger than those from the smallest release boxes, despite the 800% increase in initial area.
Sensitivity to initial area of release was found to vary throughout the NSR. Particles launched from the west of the route (which had the largest footprints) showed less variability than those in the east. The seven western release sites (Murmansk-Laptev Sea) had an average increase of 12% between the largest and smallest release sites. This contrasted with a 48% average increase in the easternmost release sites (New Siberian Islands to Bering Strait, excluding the New Siberian Islands site). In short, the larger the footprint size, the more sensitivity to the size of the initial release box.

Discussion
Understanding the advective pathways in the Arctic Ocean is an essential prerequisite for understanding the long-term fate of oil (or other pollutants) that could be released into the Arctic Ocean. Although interest in using the Arctic Ocean for commercial shipping is already increasing, particularly along the Northern Sea Route (NSR) (Liu & Kronbak, 2010;Schøyen & Bråthen, 2011), the Arctic remains a unique logistical challenge due to its harsh environment and relative lack of infrastructure and marine services (Ho, 2010). This lack of infrastructure, coupled with the often-considerable distances to ports along the NSR mean that should an oil spill occur, cleaning up an oil spill in the Arctic is more challenging than in other parts of the world. Therefore, it is necessary to consider the risk of significant amounts of oil not being recovered before the winter freezeup begins.
Should this happen, oil is likely to remain an active pollutant for a number of years due to the slowed biodegradation in the cold Arctic waters (Fingas & Hollebone, 2003). This presents a serious risk for the polar ecosystem: the food web in the Arctic is short and therefore contamination of one species can strongly affect the whole ecosystem-trophic interactions in the Arctic are simple compared to in other parts of the world, so population changes in just one key species would have cascading effects throughout the whole ecosystem (Hop & Gjøsaeter, 2013;Palumbi et al., 2008). Previous research has worked toward a quantitative risk assessment of Arctic oil spills, and the locations reached by spilled oil have been identified as key variables for understanding the ecological impacts of a potential spill (Nevalainen et al., 2017). This means that understanding the circulation patterns in the Arctic is of key importance, and it makes the advective footprints technique used here a particularly powerful tool for assessing the long-term impact of potential Arctic oil spills. Note. The metrics described in the previous 3 section (area of footprint after 9 months, straight-line distance traveled in 9 months, and fraction of particles below 100 m deep after 9 months) are recalculated with respect to year and season of release.
When considering the long-term fate of spilled oil, we need first to know how far it is likely to be transported during some specified time-frame. We took 9 months as the likely time between a spill occurring and it becoming accessible after the next melt season, and 2 years to provide a longer-term outlook, and then measured how far particles were from their start points at these intervals. Within 9 months, particles traveled a mean distance of 1,497 km from their start point (calculated along the full path). The straight-line distance showed considerable variation. A 628 km was the average for all release sites, but there was a clear east-west divide: particles from the west of the route (Barents, Kara, and Laptev) seas traveled almost twice as far as those in the east (East-Siberian and Chukchi Seas, excluding the Bering Strait): 771 km compared to 382 km. This suggests that the location of a spill is key for determining how far oil could be transported. After 2 years, particles were on average 1,223 km from their start point (484-1,662 km, depending on release location). Although the circulation in the Arctic is relatively slow compared to the rest of the world ocean, this demonstrates the first key result of this experiment: pan-Arctic modeling is required to understand the long-term (order of years) fate of spilled oil.
Having quantified the distance that spilled oil could travel, we move on to addressing the question of where it will go to and the uncertainty in how it will spread. To do this, we introduced the concept of advective footprints. This technique has recently been used in other parts of the world to study changing circulation patterns and potential impacts for environmentally sensitive areas (Polyakov et al., 2017; van Gennip et al., 2017). The trajectories from our Lagrangian experiments are in agreement with circulation pathways described in the previous literature (e.g., Aksenov et al., 2011;Carmack & Wassmann, 2006;Williams & Carmack, 2015). Known pathways are clearly reproduced, with key features such as the Arctic Circumpolar Boundary Current (ACBC) (Aksenov et al., 2011) clearly visible in trajectories launched from all along the Northern Sea Route. These, along with the Transpolar Drift Stream (Steele et al., 2004) and Pacific inflow  are the main large-scale features seen in our experiments. More detailed descriptions of specific pathways are presented in section 3.2.1.
We defined the metric ''area of advective footprints'' to quantify the uncertainty associated with these pathways. We divided the Arctic into a regular grid and counted the number of grid cells that were occupied by particles a set time after particles were released (full details in section 2.5). The sum of the areas of occupied grid cells was taken to be the area of the advective footprint. In the context of an oil spill, the size of a footprint from an individual experiment represents the spreading of the spill, while the area of the ''envelope'' of all footprints corresponds to the uncertainty in where it will go-and hence the area of ocean that would need to be searched to find it. Location of oil has been identified as a key variable for environmental risk assessments (Nevalainen et al., 2017), so it is necessary to quantify the uncertainty associated with where oil is likely to end up. Advective footprints provide a robust mechanism for doing that.
It was found that after 9 months, a typical footprint had an area of 45,762 km 2 , with a standard deviation of 13,307 km 2 , considering all release sites with an ensemble of particles initially distributed over an area of 100 km 2 . In terms of an oil spill, this figure corresponds to the area of ocean that would likely need to be affected in the spring following an unsuccessful clean-up from the previous year. This also highlights the areas likely to be at risk of potentially acute biological impact.
The size of footprints was found to be sensitive to release location, and using the results presented in section 3.2.2, we can define three main areas along the NSR, each with its own characteristic behavior: a western region (in the Barents and Kara Seas, release sites 1-7), a central region (Laptev Sea and the East-Siberian Sea, release sites [8][9][10][11], and an eastern region (in the Chukchi Sea and Bering Strait, release sites [12][13][14][15]. The west of the NSR was associated with the largest footprints after 9 months-up to 78,608 km 2 (from the Kara Sea) and on average 61,361 km 2 . The middle of the route was associated with the smallest footprints, on average 22,739 km 2 . The eastern group of release sites had typical footprints of 41,541 km 2 after 9 months. In addition to considering the spread from individual experiments, the area covered by all 30,000 trajectories from a given site (all years/seasons of release) was calculated. This metric, termed the area of the ''envelope'' of footprints represents the uncertainty in where the oil will go-i.e., the area of ocean that would need to be monitored the next spring in order to find the spilled oil. This demonstrates another important conclusion: the uncertainty in where oil could go is highly dependent on where a spill occurs. Particles from the west of the NSR tend to follow the Arctic Circumpolar Boundary Current (ACBC). They travel the farthest and are associated with the largest envelopes of footprints, and therefore a spill here would likely to require a large area (over 800,000 km 2 ) of ocean to be searched. Particles from the Siberian and eastern groups of release sites did not travel as far, and were associated with smaller envelopes of footprints-especially in the case of the East-Siberian Sea, where the envelope was only 317,615 km 2 -less than half the size of that of the largest envelopes. In short, the impact of advection on an oil spill occurring in the east of the NSR is likely to be more predictable than for a spill occurring in the west.
In addition to horizontal advection, it is important to consider the risk of potential subduction of oil. After the Deepwater Horizons spill in the Gulf of Mexico, unrecovered oil became suspended in deep water layers before eventually settling on the ocean floor (Valentine et al., 2014). We considered potential surface spills from shipping accidents in this research, but if some of this oil becomes entrained into deeper water, it is also likely to be irrecoverable. Therefore, it is important to understand where subduction occurs. We considered the fraction of particles that had sunk below 100 m depth after 9 months and also after 2 years. This was found to be of greatest concern for particles launched from the west of the NSR, particularly for release sites in the Barents Sea (see section 3.2.4).
From the three release sites in the Barents Sea (Murmansk, Barents, and Novaya Zemlya), 16% of particles were below 100 m depth after 9 months, rising to 28% after 2 years. A negligible (<1%) fraction of particles from the central group of release sites (in the Laptev and East-Siberian Seas) experienced subduction. The eastern group showed only modest subduction, mostly from particles released in the Bering Strait: 6% of particles from this site were below 100 m after 9 months and 10% after 2 years. Subduction primarily occurred as downwelling at the Eurasian shelf slope. This is consistent with results from previous investigation into downwelling in the Arctic Ocean. (Shapiro et al., 2003).

Limitations and Future Work
We are limited by dependency on a single simulation with the NEMO model: though this is a leading edge high-resolution model, intercomparison work would provide added verification and estimates of model structural uncertainty. Limitations associated with this configuration of NEMO, such as its lack of tidal forcing, means that areas strongly affected by tides may be less well represented than the Arctic in general (Padman & Erofeeva, 2004). Furthermore, despite its 1/128 resolution, given the small internal Rossby radius of deformation in the Arctic (Nurser & Bacon, 2014), NEMO is only eddy-permitting rather than fully eddyresolving in all of the Arctic. As ever, this means that higher resolution would resolve more physical processes and thus potentially provide better results.
This study is also limited in that we are only considering the portion of the oil which remains mixed into the water column and does not become trapped in sea-ice-as can happen when oil becomes encapsulated into growing ice (Afenyo et al., 2015). This can, in principle, be addressed with a similar study where the Lagrangian software used here is modified to track ice drift rather than ocean currents. However, realizing this would require significant modification of the particle-tracking software used here that is beyond the scope of this research. Additionally, oil would only be encapsulated in ice during the winter, and deciding when to follow ice or waters is itself a nontrivial question. This presents an interesting opportunity for further research, in which similar analysis using the concept of advective footprints could be put to use.
As discussed in section 3.2.6, interesting small-scale changes in current direction were seen from certain release sites in these experiments. Some can be attributed to processes described in previous literature  while others, such as a potential wind-driven blocking in the Kara Gate, provide scope for further investigation.
While this study has focused upon considering the long-term fate of spilled oil, this is not the only pollutant that is of potential interest. In principle, a similar approach could be used to assess the impact of advection on other Arctic pollutants. As with the encapsulated oil example, this may require modification of the Lagrangian software to realize this. The Arctic Monitoring and Assessment Programme (AMAP) identified four main categories of pollutant which pose a threat to the Arctic: persistent organic pollutants (POPs), ''chemicals of emerging concern'' (including flame retardants and pesticides), heavy metals and radioactive waste (AMAP, 2015). Some of these (e.g., riverine pollutants and pollution from nuclear submarines) are well suited to Lagrangian modeling as they spread from an easily defined source. This approach is not applicable in cases where we cannot accurately predict the locations where they interact with the ocean.

Conclusions
Circulation pathways and associated time scales in the Arctic Ocean have been investigated, and impact of oceanic advection for potential oil spills from the Northern Sea Route (NSR) has been explored. It has been demonstrated that pan-Arctic consideration of an oil spill's fate becomes necessary after a relatively short (1-2 years) time scale. The circulation patterns simulated are in agreement with observations.
Three main regimes have been defined to describe different sections of the NSR: western, central, and eastern. These are predominantly controlled by Atlantic Water, interior shelf-sea dynamics, and Pacific Water, respectively. Spills occurring at the west of the NSR are likely to travel the furthest (on average, 771 km within 9 months) and have the largest uncertainties associated with their pathways. Spills from the eastern and central groups travel less far (382 km), though those in the east have more uncertainty associated with their pathways than those in the center.
''Advective footprints'' were introduced to quantify the area of ocean that would likely be affected in spring following an unsuccessfully cleared up oil spill the previous season. On average, this figure is 45,762 km 2 , but it is notably higher in the western section (61,361 km 2 ) and notably lower in the central section (22,739 km 2 ) of the NSR. The 'envelope' of these footprints was introduced to quantify the uncertainty in where the oil could go-and hence how much of the ocean would need to be searched/considered potentially at risk. For a typical Arctic oil spill, this area was 676,917 km 2 .
Subduction of oil, potentially leading to unrecoverable pollution, was identified as a risk for oil spills in the Arctic, especially those occurring in the west of the NSR. It also poses a risk in the east, though to a lesser extent. Subduction is unlikely to be a concern in the Laptev or East-Siberian Seas.
This study has provided a broad overview of the circulation features that would be likely to play a significant role in the event of an oil spill from the Northern Sea Route. We have been able to provide a general description of the directions that spilled oil would likely to be carried in, and the uncertainties associated with different regions. This provides the groundwork for more focused case studies in the event of an actual spill occurring.