Typical and anomalous pathways of surface-floating material in the Northern North Atlantic and Arctic Ocean

Surface waters of the oceans carry large amounts of material, including sediment grains, plankton organisms, and ice crystals, as well as pollutants, e.g., oil and plastic. Transport and spatio-temporal distribution of this material depend on its properties and on the dynamical processes in the ocean mixed layer—currents, waves, turbulence, and convective mixing—acting at a wide range of scales. Due to its importance for marine physics, biogeochemistry and ecology, substantial research efforts have been invested in recent years in observations and modelling of ocean material transport, especially in the context of marine plastic pollution. Nevertheless, many important questions remain unanswered. In this work, numerically simulated trajectories of surface-floating particles in the period 1993–2020 are used to analyse typical and anomalous transport pathways in the northern North Atlantic and the Arctic Ocean. Model validation is performed based on additional simulations of 387 buoy tracks from the International Arctic Buoy Programme in the years 2014–2020. The trajectories are computed based on surface currents from a hydrodynamic model and Stokes drift from a spectral wave model. It is shown that due to high amplitudes of Stokes drift (comparable with wind-induced currents in ice-free parts of the domain of study), combined with high directional variability, the drifting paths are substantially modified in ice-free regions, underlying the important role of wave-induced currents in surface material transport. A statistical analysis of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 1.6\, 10^8$$\end{document}∼1.6108 trajectories reveals patterns of connections between nearshore locations in the domain of study, the associated drift times and path sinuosity. Seasonal variability of transport, which differs between the Arctic Ocean and the North Atlantic, is found for typical transport routes following the larger-scale circulation patterns. Crucially, in both sub-domains episodic, but very strong transport events between otherwise isolated locations occur, associated with anomalous atmospheric circulation and, arguably, providing ‘windows of opportunity’ for dispersal of various organisms to new locations. It is shown for two examples in the North Atlantic region that an unusual combination of atmospheric circulation indices explains the anomalous transport, thus providing a predictive tool for future events. In the Arctic, analogous phenomena are modified by the state of the sea ice cover.

Supplementary Note S1. The role of the diffusion coefficient K h As noted in the main text (see Methods), sensitivity simulations have been performed for trajectories originating at ICE in order to assess the influence of the value of the diffusion coefficient K h on the modelled particle trajectories. Three values of K h have been tested, K h = 29 m 2 ·s −1 (used in the simulations described in the main text), K h = 290 m 2 ·s −1 and K h = 2.9 m 2 ·s −1 . Supplementary Fig. 3 shows the corresponding probabilities of drifters reaching from ICE to the remaining 15 destination regions (see also connection matrix in Fig. 3a in the main text). As can be seen, for all connections with probabilities higher than 0.1%, their values obtained with the three different K h s are very similar and remain at the level below 1%. Only connections to FJL, GLN and SBW show larger differences, but those destinations are reached by too few drifters for the statistics to be robust (this is the reason those connections are not analyzed in the discussion in the main text). Analogously, the statistics of the minimum and median drift times for connections from ICE remain largely unaffected Supplementary Fig. 4), with the same exceptions as those described above -although, as expected, overall larger K h lead to slightly shorter times, and the differences in times obtained with different K h get more pronounced for longer trajectories, e.g., those to BEI or SVA.
The time series of release and landing events, shown for two examples (NOR and SVA) in Supplementary  Fig. 5, exhibit some differences for different K h -again, as expected, the peaks in those time series are generally smaller for smaller K h and vice versa -but, crucially, the differences are very small between K h = 29 m 2 ·s −1 and K h = 290 m 2 ·s −1 , the two values that can be seen as limits of K h reported in the literature.
Overall, the results of the sensitivity tests show rather limited influence of K h on the statistical properties of the modelled trajectories, suggesting that the results obtained with K h = 29 m 2 ·s −1 , discussed in the main body of this paper, can be treated as representative for a wide range of realistic values of the diffusion coefficient. Supplementary

Supplementary Note S2. Validation against buoy data
As mentioned in the main text, buoy tracks from the International Arctic Buoy Programme (https://iabp. apl.uw.edu/) are used here as a source of in situ data for validation of the model. The usage of IABP data has several limitations due to the fact that the buoys' drift speeds and directions are a net result of forces from the uppermost ocean and the lowermost atmosphere -they do not strictly fulfill our assumptions of drifters subject to forcing at depth z = 0. Moreover, the influence of the wind and current forcing on the buoys tracks is very hard to estimate, as several types of buoys are used, with different mass, shape, and oceanographic/atmospheric instruments attached to them. Also, the role of winds and currents certainly changes with changing ice conditions (ice type, concentration, etc.). Nevertheless, the IABP buoy tracks are routinely used in analyses of sea ice drift patterns, dispersion etc., not only in the central Arctic, but also in very dynamic regions as e.g. the Fram Strait/East Greenland Current. They are also used for validating satellite products of ice motion (including the data assimilated in HYCOM). Therefore, keeping the above limitations in mind, we are convinced that the comparison of our trajectories with the IABP data provides useful insights into the quality and reliability of our results.
In the following analysis we use LEVEL 2 data (i.e., after basic quality control), available in the years 2015-2020, supplemented with LEVEL 1 data from 2014 in order to increase the number of data points (and, in particular, the number of tracks in the open-water North Atlantic part of the domain of study). The time series of buoys' positions are interpolated in time onto the 12:00 hour of each day, and those daily time series are used for model validation. After removing tracks with unrealistically high displacements and/or obviously incorrect positions, a total of 387 tracks are retained for further analysis, ranging from 1 to 39 months in duration (with a median value of 209 days) and from 1 to 15000 km in length (with a median value of 1960 km). The initial and final positions of those tracks are shown in Supplementary Fig. 6. For each track, the model is run 20 times (10 times without the Stokes drift, and 10 times with Stokes drift taken into account), initialized from a location randomly selected from within a square region with side length equal to 5 km, centered at the initial buoy location (i.e., according to the same procedure as used in the main simulations analyzed in the text of the paper). As the main measure of the model performance, the relative error d best /L has been selected, where d best denotes the distance (in meters) between the simulated and observed final positions of the buoy, computed for the 'best' path, and L denotes the total length (in meters) of the observed buoy path. The ratios d best /L are computed separately for simulations with and without Stokes drift. The summary of the d best /L statistics is shown in Supplementary Fig. 7 in the form of the respective pdfs and cdfs. As can be seen, in ∼50% of the cases the simulated final positions have errors smaller than 10%, and for close to 80% of cases the error is smaller than 20%. In general, the pdfs of d best /L can be well described by exponential curves. Notably, when all trajectories are taken into account, the role of the Stokes drift vanishes -the exponential fits for the results obtained with/without the Stokes drift are essentially the same. The reason for that is the fact that the majority of IABP buoys are located in sea ice, when the role of waves, and thus also the role of the associated Stokes drift, is negligible. When analogous pdfs/cdfs are computed for open-water cases (i.e., trajectories located in the North Atlantic or in the ice-free regions of the Arctic Ocean during the summer), the model performance is significantly better when the Stokes drift is taken into account (compare the blue and red curves in Supplementary Fig. 7c.d).

Supplementary
Importantly, as well, the bias of the length of simulated trajectories for simulations with Stokes drift is very small: the mean normalized difference between the total modelled and observed path length, (L s −L)/L, equals −0.07, and the corresponding standard deviation 0.56. The bias for model version without Stokes drift is much higher and equals −0.26, in accordance with the fact that the paths without the contribution of wave-induced drift are usually shorter. The scatter in that case is similar and equals 0.53.
In order to get a better insight into the model performance and behaviour, we show selected examples of trajectories in Supplementary Figs 8-14. In all figures, the initial position and the trajectory of the IABP buoy are shown with a yellow triangle and a red curve, respectively; the simulated trajectories and their final positions are marked with magenta curves and white circles (version without Stokes drift) and with blue curves and blue circles (version with Stokes drift), respectively. The observations made based on those examples (and several other tracks, which are not shown) are as follows: In most parts of the central Arctic Ocean, the model accurately reproduces both drift speed and direction of the buoys ( Supplementary Fig. 8), presumably thanks to the assimilation of satellite sea ice drift data used in the ocean circulation model, on which our simulations are based.
In the Chukchi and East Siberian Seas, and in particular in the region around the Wrangel Island, the trajectories often have very irregular shapes, with several loops and rapid shifts of direction, combined with a relatively low average drift speeds. Although the modelled trajectories generally reproduce those overall path characteristics, there is a large scatter between individual paths, especially during periods of weak atmospheric/oceanic forcing and low sea ice concentration. Nevertheless, in most cases the mean drift direction is correctly reproduced (Supplementary Fig. 9).
The simulated trajectories along the east coast of Greenland usually reveal large scatter in terms of their landing locations, as the 'fate' of individual paths is very sensitive to their distance from shore and small perturbations in forcing. However, in most cases (including the examples in Supplementary  Fig. 10), some of those trajectories reach land very close to the observed landing location. In this case, as in the one described in the previous point, the drift speeds computed without Stokes drift generally tend to be too slow compared to observations. The model captures also drift paths within the cyclonic recirculation region to the south of the Denmark Strait ( Supplementary Fig. 11).
Reproducing trajectories during weak and/or highly variable wind and wave forcing is a challenge ( Supplementary Fig. 12 shows an example from winter-time in the Barents Sea). However, by definition, those situations contribute very little to the long-distance transport between coastal regions, which is the main subject of the paper. and for buoy tracks in ice-free areas (c,d): probability density functions (a,c) and cumulative probability distributions (b,d) of the relative model error, defined as the ratio between d best (defined as the minimum distance, among all simulated paths, between the final predicted and observed buoy position) and the total path length L. Note that both versions of the exponential fit in (a,b) are essentially identical, so that only one fit is shown. Bin width and number of data points in (a,b) equal 0. 05 and 387, respectively; in (c,d), the values are 0.1 and 56. See text for more details.
In regions without sea ice and/or with low ice concentration, the Stokes drift enhances the average drift speeds ( Supplementary Fig. 13a,b) and contributes to the landing of particles, as in the example in Supplementary Fig. 13c, where none of the no-Stokes trajectories, but nine out of ten of the Stokestrajectories reached land (some in vicinity of the observed landing location).
Finally, there are several cases when the model failed to reproduce the observed trajectories, as in the examples in Supplementary Fig. 14. Those very large discrepancies are typically associated with model inaccuracies in regions located at the boundaries of two neighboring circulation cells, like e.g., the Beaufort Gyre and the Transpolar Drift leading to the Fram Strait ( Supplementary Fig. 14b). As the statistics presented earlier in this section shows, those cases constitute only a small amount of the total number of trajectories.