1 Introduction
The Lagrangian study of transport and mixing in the ocean is of fundamental interest to ocean modellers (van Sebille et al., 2018, 2009; LaCasce, 2008). In particular, the analysis of data obtained from Lagrangian drifting objects greatly contribute to our knowledge of ocean circulation, e.g. through analysing the accuracy of numerical and stochastic models (Huntley et al., 2011; Sykulski et al., 2016), or the use of the data to better understand various pathways and where to search for marine debris (Miron et al., 2019; van Sebille et al., 2012; McAdam and van Sebille, 2018).
Meehl (1982) used shipdrift data to create a surface velocity data set on a grid. These velocities were used to simulate the Lagrangian drift of floating objects in Wakata and Sugimori (1990). More recent works focus on using drifting buoys to derive Lagrangian models to discover areas where floating debris tends to end up (van Sebille, 2014; van Sebille et al., 2012; Maximenko et al., 2012). Advances in technology have resulted in much better data quality, which now permits the use of more detailed methodology. The newer models provide densities of where debris ends up on grid scales as small as .
In this paper, we propose a novel computationally fast method for estimating a so called “most likely pathway” between two points in the ocean, alongside an estimated travel time for this pathway. The method is purely datadriven. We demonstrate our methodology on data from the Global Drifter Program (GDP), but the method is designed to work with any drifter data set. Additionally, we develop and test related methodology for providing uncertainty on both the pathways and the travel times. Our method is automated with little expert knowledge needed from the practitioner. We provide a set of default parameters which allow the method to run as intended. The method simply takes in a set of locations within the ocean, and outputs a data structure containing most likely paths and corresponding travel time estimates between all pairs of locations. We focus on a global scale: we aim to provide a measure of Lagrangian connectivity for locations which are thousands of kilometres apart. Individual drifters are unlikely to connect two arbitrary locations far apart, hence the need for our methodology which fuses information across many drifters.
A tool which predicts travel times is of practical use in many fields. For example in ecological studies of marine species, genetic measurements are taken at various locations in the ocean. Euclidean distance is often used as a measure of separability and isolationbydistance (Becking et al., 2006; Ellingsen and Gray, 2002) to find correlations with diversity metrics or genetic differentiation between communities or populations of organisms. Due to various currents and land barriers, we expect Euclidean distance to often be a poor measure of how ‘distant’ or dissimilar two points are. The method proposed in this work would use the estimated travel times to supply a matrix containing a Lagrangian distance measure between all pairs of locations. This matrix can then be contrasted with a pairwise genetic distance matrix between these locations and will yield new insights. In many instances the Lagrangian distance matrix will be more correlated with genetic relatedness than a Euclidean distance matrix. This observation was already made in the Mediterranean Sea when studying plankton (Berline et al., 2014), and off the coast of California for a species of sea snail (White et al., 2010). Both of the works by Berline et al. (2014) and White et al. (2010) rely on simulating trajectories from detailed ocean currents data sets to estimate the Lagrangian distance. This methodology does not scale globally and relies on simulated trajectories from currents rather than real observations.
In Figure 1, we show seven locations plotted on a map with ocean currents. We shall use these locations as a proofofconcept example throughout this paper. The exact coordinates are given in Table 1 in Section 5. The aim is to introduce and motivate a method which provides an estimate as to how long it would take to drift between any two of these locations. For example, the travel time from location 2 to location 3 in the South Atlantic Ocean should be smaller than the return journey due to the Brazil current. We choose to include locations in both the North and South Atlantic as we wish to demonstrate that the method successfully finds pathways linking points which are extremely far apart.
1.1 Comparison with Related Works
In this section we give a brief overview of techniques that have used the Global Drifter Program to achieve a similar or related task. The work by Rypina et al. (2017)
proposes a statistical approach for obtaining travel times. A source area is defined such that at least 100 drifters pass through the source area. The method focuses on obtaining a spatial probability map and a mean travel time for every
bin outside of the source area. This method successfully combines many trajectories, however for multiple locations one would have to decide on a varying grid box for each location of interest. Such a grid box must be manually chosen by the practitioner meaning that the method does not scale well with multiple locations. Rypina et al. (2017) also focus on estimating a mean travel time, where our method provides a travel time associated with the most likely path, and is hence more akin to estimating a mode or median travel time.The method by van Sebille et al. (2011), which proposes the use of Monte Carlo Super Trajectories (MCST), could naturally be used to estimate travel times. This method simulates new trajectories as sequences of unique grid indices along with corresponding travel time estimates for each part of that journey. The method is purely data driven i.e. they only use real trajectories to fit the model. The travel time and pathway we supply here should be similar to the most likely MCST to occur between the two points. The advantage of our methodology is that we do not base the analysis on a simulation, such that the results from the method described in Section 3 are not subject to any randomness due to simulation.
Various other works have made attempts to compute Lagrangian based distances. For example, Berline et al. (2014) used numerically simulated trajectories to estimate Mean Connection Times within the Mediterranean Ocean. Smith et al. (2018) used MCST to estimate various statistics of how seagrass fragments could drift from the South East coast of Australia to Chile. Specifically, Smith et al. (2018) simulated 10 million MCST starting from the SE coast of Australia and only 264 (0.00264%) of the simulated trajectories were found to travel roughly to the Chilean coast.
The approach by Jönsson and Watson (2016) uses simulated drifter data to construct connectivity matrices between locations in the ocean. As the matrix is sparse, Dijkstra’s algorithm is used to connect arbitrarily distant locations in the ocean to measure Lagrangian distance. Although this method may at first glance bare similarities with our method (specifically in the use of Dijkstra’s algorithm), there are in fact many differences. First of all, the method uses simulated trajectories whereas we use real drifter trajectories. Secondly, Dijkstra’s algorithm is performed by Jönsson and Watson (2016) on the connectivity matrix (which finds minimum connection times between locations), whereas our approach uses Dijkstra’s algorithm on the transition matrix which describes a probabilistic framework for drifter movement. We found the latter approach to perform much better with real data. Finally, we cannot directly implement the approach described in Jönsson and Watson (2016) as only connectivity values higher than one year are used. For real data such a step would result in a very sparse connectivity matrix making the method infeasible. An initial analysis we conducted using similar methodology achieved poor results.
In contrast to all these previous works our methodology provides three novel contributions: (1) a computationally efficient approach for simultaneously finding most likely paths and travel times across multiple locations without requiring simulated trajectories; (2) the use of the (H3) spatial indexing system for discretization of drifter data; and (3) methods to address error due to grid discretization and the uncertainty from sparsity of observations.
The method we propose is computationally efficient even with a large set of locations. In contrast, if we used MCST as in Smith et al. (2018), 10 million trajectories would be released from each location of interest to obtain an estimate of travel time to all other locations. This procedure would be required for each location of interest resulting in a very large number of trajectories being computed. In the method we propose we only need one run per location of an efficient shortest path algorithm which may run in a matter of seconds. Also, as we do not rely on simulated data: if it is found that an area is not accessible by our method (i.e. there exists no pathway), that means that there is insufficient data in the drifter data set to access that point. However, in a simulation approach, the pathway may not have been generated across the simulations, even though there was in fact sufficient evidence in the data for one to exist, resulting in potentially missed pathways.
1.2 Structure of the Paper
The structure of the paper is as follows. In Section 2, we give a brief overview of the Global Drifter Program and the movement of the objects which the program tracks. In Section 3, we introduce the methodology which obtains the most likely path and a corresponding travel time. We provide guidelines on how to quantify uncertainty due to randomness and disretization errors in Section 4. Results using the locations of Table 1 are given in Section 5. Finally, we summarise our findings in Section 6.
2 Background and Notation
2.1 Global Drifter Program
The Global Drifter Program(GDP) is a database managed by the National Oceanographic and Atmospheric Administraction (NOAA) (https://www.aoml.noaa.gov/phod/gdp/) (Lumpkin and Pazos, 2007). This data set contains over 20,000 freefloating buoys temporally spanning from February 15, 1979 through to the current day. These buoys are referred to as drifters. In this analysis we use 24,030 drifters; the spatial distribution of which are shown in Figure 2. The drifter design comprises of a subsurface float and a drogue sock. Often this drogue sock detaches. We refer to the drifters which have lost their drogue as nondrogued drifters, and drogued for those which still have the drogue attached.
Here we use the drifter data recorded up to August, 2019. We use data which has been recorded from drogued drifters only. This is not a restriction, as it would be straightforward to simply use the data from nondrogued drifters if a practitioner was interested in a species or object which experiences a high wind forcing, or a combination of both if it is believed that the species followed a mixture of near surface and windforced currents. The data is quality controlled and interpolated to six hourly intervals using the methodology from Hansen and Poulain (1996). These interpolated values do contain some noise due to both satellite error and interpolation, however, the magnitude of this noise is negligible in comparison to the size of grid we use in Section 3. Hence, we ignore this noise and treat the interpolated values as observations. For the same reason we note that the interpolation method used is not important here, instead of the six hourly product we could use the hourly product produced by methodology proposed by Elipot et al. (2016).
The value of using the Global Drifter Program is we obtain a true modelfree representation of the ocean. All phenomena which act on the drifters are accounted for in the data set. The other common approach is to first obtain an estimate of the underlying velocity field, then simulate thousands of trajectories using the velocity field. While this simulation approach is often satisfactory in some applications, the models generally do not agree completely with the actual observations.
2.2 Notation
Here we use to be a geographic coordinate corresponding to latitude and longitude respectively. We refer to the longitudelatitude grid system using the notation , which means each grid box goes along the longitude axis and along the latitude axis. We use bold font for any data which is in longitudelatitude pairs; i.e , and nonbold text for either a grid index or a single number. We use to denote the set of all possible grid indices.
2.3 Capturing Drifter Motion
We define the drifter’s probability density function as
where the drifter started at at time and moved to position at time , where and are longitudelatitude pairs. In the absence of a model, this probability density cannot be estimated continuously from data alone. Therefore, we follow previous works which discretize the problem (Maximenko et al., 2012; van Sebille et al., 2011; Miron et al., 2019; Rypina et al., 2017). Instead of considering , we consider where is some set of states which correspond to a polygon in space; we will define how these are formed in Section 3.2. Often these states are simply degree boxes (e.g. as used in Figure 2). As in Maximenko et al. (2012), we assume that the process driving the drifter’s movement is temporally stationary. That is:
i.e. the probability of going from to depends only on the time increment. The probability does not depend on the start or finish time.
Furthermore, given that we are using data which is observed at regular and discrete times, we shall only consider discrete values of time. Let be a sequence of random locations which are interested in, each entry can take the value of anything within . We define the probability as the probability that the position at time is given that the state at time was where .
A Lagrangian decorrelation time causes the drifter to ‘forget’ its history (LaCasce, 2008). We aim to chose a quantity which is globally higher than the Lagrangian decorrelation time. We call this quantity the Lagrangian cut off time . The reasoning behind using this time is that if we consider a sequence of observations at least apart then the following Markov property is satisfied:
(1) 
The first term in Equation 1 is the probability that both is and is . In other words, the Markov property states that the probability of transitioning to a state at time is independent of all the past states at times and earlier, given the state at time is known. In this case, the physical time difference associated with and being larger than the chosen Lagrangian time scale , validates the use of the Markov assumption.
For the rest of this paper we assume that the time between observations is at least . Which allows us to use the Markov property from Equation 1
freely. In so doing, alongside the simplification of discretizing locations, this allows the problem to be treated as a discrete time Markov chain. Here we fix
days as this matches previous similar works (Maximenko et al., 2012; Miron et al., 2019). The estimated decorrelation time for the majority of the surface of the Ocean is likely to be lower than 5 days (e.g. see Zhurbas and Oh (2004) for the Pacific and Lumpkin et al. (2002) for regions in the Atlantic).3 Method for Computing the Most Likely Path and Travel Time
Maximenko et al. (2012) and van Sebille et al. (2012) focus on the use of a transition matrix estimated from drifters to discover points where drifters are likely to end up. In this section we build on such an approach by providing a method to take such a matrix and provide an ocean pathway and travel time.
In Section 3.1, we explain in detail how the transition matrix is formed. As a grid system is needed to form the discretization of data we introduce our chosen system in Section 3.2. Then in Section 3.3, we describe how we estimate the most likely path of a drifter to have taken. Finally, in Section 3.4 we explain how we turn the most likely path and transition matrix into an estimate of travel time. We give a summary of how this all comes together in the pseudocode in Algorithm 1.
3.1 Transition Matrix
The location of a drifter at any given time is a continuous vector in
, the longitude and latitude of the point. We define an injective map which maps this continuous process onto a discrete set of states which are indexed by integers in . We define the map as follows:(2) 
We aim to make a Markov transition matrix of size rows and columns. Where denotes, the probability of moving from to in one time step. Very similar to the approach in Maximenko et al. (2012), we form our transition matrix using a gap method. In each drifter trajectory we only consider observations as a pair of points days apart. When using this method for other applications we advise using to be higher than the decorrelation time of velocity to justify the Markov assumption.
Consider a trajectory as a sequence of positions where is the out of trajectories, is the number of location observations in the trajectory, and are the longitudelatitude positions. First, we map each trajectory into observed discrete states. We will denote these states as follows,
For each we estimate the relevant entry of our transition matrix through using the following empirical estimate:
(3) 
Note that we take gaps of as observations are every 6 hours in the GDP application. We expect that states in which are not spatially close will have nonzero entries, therefore, the matrix will be very sparse.
3.2 Spatial Indexing
Clearly the resulting transition matrix described in Section 3.1 strongly depends on the choice of grid function in Equation 2. Most previous works (McAdam and van Sebille, 2018; van Sebille et al., 2012; Rypina et al., 2017; Maximenko et al., 2012) use longitudelatitude based square grids where all grid boxes typically vary between and . A grid cell around the equatorial region (0 degree latitude) will be approximately equal area to a square box. However, if we take such a grid above latitude i.e the Norwegian sea, the grid cell will be approximately equal area to a square box.
There are a few other choices which we argue are more suitable for tracking moving data on the surface of the Earth. Typically three types of grids exist for tessellating the globe: triangles, squares, or a mixture of hexagons and pentagons. Here we choose to use hexagons and pentagons as they have the desirable property that every neighbouring shape shares precisely two vertices and an edge. This is different to say a square grid where only sidebyside neighbours share two vertices and an edge, whereas diagonal neighbours share only a vertex. This equivalence of neighbors property for hexagons and pentagons is clearly desirable for the tracking of objects as this will result in a smoother transition matrix.
We specifically use the grid system called H3 by UBER (UBER, 2019). This system divides the globe such that any longitude and latitude coordinate is mapped to a unique hexagon or pentagon. This shape will have a unique geohash which we can use to keep track of grid index. The benefit of using such a spatial indexing system is that attention is paid to ensuring that each hexagon is approximately equal area. We use the resolution 3 index where each hexagon has an average area of . A square box of size has roughly the same area as this which is very similar to the size of a grid cell near the equator. An example of an area tessellated by H3 is shown in Figure 3. Other potential systems which could be used include S2 by Google which is a square system, or simply using a longitudelatitude system as various other works do.
3.3 Most Likely Path
For our analysis, the first step is to define a most likely path. A path is simply a sequence of states such that the first element is the origin and the last element is the destination. We also require that two neighboring states are not equal to each other.
Definition 1 (Path).
We define the space of possible paths as between the origin and destination as the following:
(4) 
With a cardinality operator which denotes the length of the path.
Given the transition matrix we define the probability of such a path:
(5) 
Definition 2 (Most likely path).
Consider any path . By the most likely path we mean the path which maximises the probability of observing that path.
(6) 
Optimising Equation 6 appears intractable at first glance in its current form. However, we consider the following form for :
Then we use the fact that:
(7) 
Now consider a network with nodes being every state in and edges being connections between these nodes. Two nodes, say , are connected if we have that , which then give this edge a weight of . With this network setup the optimization of Equation 7 can be solved efficiently using the vast literature on shortest path algorithms (Gallo and Pallottino, 1988; Dijkstra, 1959). In our methodology we employ Dijkstra’s algorithm to solve Equation 7.
3.4 Obtaining a travel time estimate
The most likely path is often a quantity of interest in itself, however we can also obtain a travel time estimate of this path. The method should be fast and efficient as it should be able to run for large sets of locations quickly. We achieve this by giving a formula to estimate the travel time based on the transition matrix.
Consider the path which we aim to estimate the time of to be . To start we assume that if the current state is then the next state is sampled from a categorical distribution, where the parameters are simply defined by the row . The categorical distribution with parameters , , simply says the probability of drawing is .
Now we assume that the only possibility is that the drifter follows the path we are interested in. So must be followed by . Now we use to index time and suppose
, then we are interested in the random variable
where , i.e. the time which the object moves between and . Note that the only possibility for states is that they are all , otherwise the object would not be following the path of interest. Therefore, we obtain the distribution of as follows:(8) 
Note that if we set in Equation 8 we get:
(9) 
which is the probability distribution function of a negative binomial distribution with success probability
and number of failures being one. We denote the random variable for the travel time between and as . As the negative binomial distribution corresponds to the time until a failure, we are interested in taking one time increment longer than this as we require to be the time that we move from to i.e. the time of the failure. Therefore the distribution of exactly follows . Also, note that is in units of the chosen Lagrangian cutoff time .To get the expectation of Lagrangian times we consider the sum of all the individual parts of the travel times , such that we obtain:
(10) 
where we have used that the expectation of the negative binomial is .
We could attempt to obtain a simple variance estimate for the estimate
with classical statistics. However, we would only be able to account for variability in the estimates of the entries , i.e. we would have to assume is known. In our case we are interested in the time of , which is an estimate as it depends on . Obtaining any analytical uncertainty in the estimation of the most likely path would be intractable due to the complexity of the shortest path algorithm. Therefore, we propose to address the issue of uncertainty in and due to data randomness in Section 4.2 using the nonparametric bootstrap. To finish this section, we provide the pseudocode for our approach in Algorithm 1.4 Stability and Uncertainty
In this section we provide further methodology to address practical issues that may arise when applying the methodology of Section 3 to a given data set. A key issue is the final results of the algorithmic approach may strongly rely on the precise grid system chosen in Equation 2, especially in regions with irregular land boundaries. We propose a solution to the problem in Section 4.1 which is based on performing rotations of the grid system to quantify discretization error. Then in Section 4.2 we present a method to estimate uncertainty in and using the bootstrap.
4.1 Rotation
The first operation we perform in Algorithm 1 is to pass the observations through the injective function to setup the data for the discretization problem. This step is incredibly important. McAdam and van Sebille (2018) give examples where a poor choice of will cause spurious connectivity which does not exist between ocean basins. An extreme example being the method implying that a drifter is able to cross the Panama land mass. Even though there is a canal, it is almost impossible for a drifter to cross that canal and has never happened to date. There is also an inherent bias introduced by this grid step. If we are interested in two different points which fall in the same spatial grid, they are estimated to have the same travel time. Here we suggest an adaptation to the grid system step to remove such biases and errors.
For example, see Figure 4, and suppose we expected the shortest path resulting from the method described in Section 3 is to travel from the topleft grid, to the middle grid, then finally into the bottom right grid, e.g. if there was a velocity field going in the positive and negative direction. Denote the travel time as as the estimate for how long it takes to get from to using our algorithm, using the superscript to identify that the black grid is the grid system which is used for . In such a case the estimate of will be zero as and are in the same grid. However, the true difference in travel times should be positive, as is closer to and has a shorter distance to travel with the flow. Consider if we were to estimate the distance using the same method except is now the red grid system. Suppose the shortest path between and is still the same, (top left, middle, then bottom right). Let the resulting travel times be denoted . We now have that is at least , i.e. that the travel time from to is higher than the travel time from to .
We expect that if we take lots of such shifted grids that some travel times may end up too low (e.g. 0), others may be too high, and some are about right. By averaging out many shifted grids we obtain a travel time around what should be expected. So in the example given in Figure 4, we would use . In addition, by averaging estimates of the distance between and we obtain a nonzero estimate, whereas in the original method just using the black grid as this estimate would always be zero.
In Figure 4, we look at a flat simple square grid—it is very simple to perturb this grid by simply shifting it. In global applications such as the one we consider here we are using complicated spatial grid systems. Rather than trying to reconfigure the grid system, instead we suggest a more universal alternative. We suggest rotating the longitudelatitude locations of all the relevant data using random rotations. Such a strategy will work for any spatial grid system as it just involves a prepossessing step of transforming all longitudelatitude coordinates. Note that for each rotation we are required to reassign the points to the grid and reestimate the transition matrix. These are the two most computationally expensive procedures of the method. For the rest of this section we briefly summarise the random rotation on a sphere algorithm proposed by Arvo (1992) which we employ.
To start, we make the simplifying assumption that the earth is spherical, and project longitudelatitude into their polar coordinates, on a unit sphere:
(11)  
With inverse
(12) 
Usually we would consider a radius term in these equations, here we consider the radius of the earth to be unit one.
In general a rotation matrix is denoted as the group SO(3), which we define as the group of matrices:
Each rotation matrix corresponds to a rotation on the sphere. There are multiple approaches to achieve this. Here we use a Euclidean angles based approach. The Euclidean angle approach is described in (Arvo, 1992); the basic idea is that first we rotate a sphere around the vertical axis, then rotate the north pole of this rotated sphere to a random position. The first step is simply a 2d rotation matrix, keeping the third coordinate fixed:
(13) 
The second rotation is not as simple. The rotation involves defining a rotation such that it will bring the north pole to any other point on the sphere with equal density. Specifically, we set our rotation matrix as follows:
(14) 
Arvo (1992) shows that if where and , that
will follow a uniform distribution over the sphere.
In summary, the steps for a random rotation involve:

Transform longitudelatitude coordinates to polar coordinates according to Equation 11.

Generate random numbers and . Use these random numbers to set according to Equation 13, and set according to Equation 14.

Set .

Set .

Obtain transformed longitude and latitude through applying Equation 12 to .
To obtain travel times which have the bias removed, we sample rotation matrices . We then run Algorithm 1, however as a prepossessing step we rotate all locations of the drifter trajectories and locations of interest. For each rotation matrix this will result in a set of distances . To obtain our actual estimate of a travel time we use a mean of all the distances:
(15) 
4.2 Bootstrap
Here we provide details of a bootstrap estimate of the uncertainty of travel times. If we required a rough estimate of uncertainty we could consider that the , the most likely path, is fixed and then estimate . However, this would be a poor estimate. The estimate assumes: (1) that the transition matrix entries follow a certain distribution, and (2) that the path is the true most likely path. In reality neither of these are true, they will both just be estimates. The transition matrix elements are estimated from limited data and the shortest path strongly depends on the estimated transition matrix, e.g. a small change in the transition matrix could result in a significantly different path. Therefore, we obtain estimates of uncertainty by bootstrapping (Efron, 1993).
Bootstrapping is a method to automate various inferential calculations by resampling. Here the main goal is estimate uncertainty around . One of the main benefits of bootstrapping is that it allows us to estimate uncertainty of the whole complicated process in Section 3. The bootstrap involves first resampling from the original drifters to obtain a new data set. We do this by resampling from the original set of drifters with replacement. This is done by sampling values with replacement from to , then for the sampled value which takes the value , set . We call a bootstrap sample. Then we use as the input to Algorithm 1.
We do this resampling times to obtain estimates of , we denote these bootstrap estimates as
. We then estimate our final bootstrapped mean and standard deviation estimates as the following:
(16) 
This method also results in sampled most likely paths . Which may be compared to the original path. Other benefits of the bootstrap include automatic highlighting of paths which rely on a small amount of data. For example, suppose the method found a direct route as the most likely path, but this path was only possible because one or two drifters filled in the transition matrix for that area. Additionally, suppose there is another more circuitous path which our method finds is less likely but many more drifters take. In many of the bootstrap estimates we will take the longer more circuitous path. Such an example may end up with a distribution of parameters with two modes. One mode at the travel time for the direct path, and another mode for the travel time for the longer path.
5 Application
Longitude  Latitude  

1  9.0  25.5 
2  25.0  5.0 
3  45.0  40.0 
4  69.0  39.0 
5  42.5  41.5 
6  42.0  27.5 
7  93.2  24.8 
We use the locations given in Table 1 for the demonstration of the method described in this paper. These locations were chosen for multiple reasons; (1) they were placed on or near Ocean currents, such as the South Atlantic current, the Equatorial current and the Gulf Stream; the magnitudes of which can be seen in Figure 1, and (2) stations were placed in both the North and South Atlantic to show how the method can find pathways which are not trivially connected. First we go over an application of the vanilla method from Section 3, then we provide brief results using the adaptations using bootstrap and rotations from Section 4 in Section 5.1 and Section 5.2 respectively.
Prior to our analysis we take two practical steps to improve the reliability of the method. The first is that we find the states corresponding to and (two points on the Panama land mass), then remove the corresponding row and column from . If this step is not taken we quite often get pathways crossing the Panama land mass, resulting in impossibly short connections to the Pacific Ocean. The second is that if for any state the denominator in Equation 3 is less than , i.e. the number of drifter observations starting in state is less than , then we drop the corresponding row and column of . The second step is not strictly necessary, but is taken to reduce spurious connections and improve the reliability of results.
As a quick note on computational speed, the method is quite fast and efficient. All using a single 2.8 GHz core on a laptop using a Python implementation of the algorithm, then the basic method described in Section 3 takes around 4 minutes to run fully and obtain all distance matrices. Around 3min 40 seconds of this is spent on passing through the spatial index and the estimation of . Estimating distances to extra locations is quite fast (less than 20 seconds for a set of new locations). Each rotation takes a full run of the original method i.e. around 4 minutes per rotation. One bootstrap sample, given that we have already matched the stations and drifters to their grid indices, takes around 30 seconds in total. The rotation methodology and bootstrap can be run in parallel easily. The transition matrix is around 200MB such that the memory usage of the method is relatively light.
Figure 5 shows the pathways between a representative sample of the stations. First we note what features are observed in the most likely path. The Gulf Stream is used on any path trying to access locations 4, 5 or 6 in Figure 5. See in Figure 5 when going from location 3 to 5 that the method chooses to enter the Gulf of Mexico and then uses the Gulf Stream to access location 5, even though the actual geodesic distance of this path is quite long. Other examples which use the Gulf Stream include and . Generally, any of the paths leaving location 1 and attempting to travel northwest uses the Benguela Current, for example Figure 5 , and .
The travel times obtained between the sample stations in Figure 5 show some interesting results regarding the lack of symmetry. When going from location 2 to location 4 we estimate quite a long most likely path in terms of physical distance. However, the resulting travel time of this path (0.5 years) is almost a quarter of the travel time from the most likely path from location 4 to location 2 (1.9 year)  which is much shorter in distance. We expect this is because the path going from location 2 to location 4 follows very strong currents such as the North Equatorial current and the Gulf Stream. Another surprising result is that going from 3 to 5 and vice versa are relatively close seeing that 3 to 5 actually uses the Gulf Stream but the return does not. In the most likely path from 3 to 5, up until around latitude the travel time is 2.9 years, which we expect as the pathway seems to be going against the Brazil current. After this point the rest of the path takes the remaining 0.8 years despite the remainder being over half the actual physical distance of the pathway. We expect this short time is due to the method finding a pathway along the North Brazil current, followed by the Caribbean current, followed by the Gulf Stream.
Figure 6 shows the travel time distribution from location 1 to the entire globe. One thing to note about this method is that the most likely path is not always the shortest path. This results in the travel time distribution not necessarily being spatially smooth. Consider the discontinuity line around 5 degrees in the Pacific ocean in Figure 6. In Figure 7 we plot the two paths, to two points which are only latitude apart, such that each one is on either side of this discontinuity. Both paths start by using the Antarctic Circumpolar Current, until we reach the middle of the Pacific ocean. We see that the path going to latitude takes a more direct approach going diagonal through the middle of the Pacific then up to . Whereas, the path going to latitude follows the South Pacific current and the Antarctic circumpolar current up to the Peru current and then the Equatorial current to reach the point of interest. The resulting path is longer in distance but significantly shorter in estimated travel time.
5.1 Bootstrap
To show the value of the bootstrap we show the results for one particular pair of stations, the pathway going from location 1 to location 3 and back. The pathways which result from the bootstrap are shown in the bottom panel of Figure 8. The darker lines on the figure imply that that this transition is used more often. We see that for most of the journey the darker lines closely follow the original path. The bootstrap discovers some slightly different paths, for example around Longitude the path going from 3 to 1 occasionally seems to find that going further south is a more likely path. Also, around the beginning of the path going from 1 to 3, we see that the most likely path taken most frequently by the bootstrap samples often does not follow the most likely path from the full data.
The main goal of the bootstrap is that we obtain an estimate of the standard errors. In this case we get a standard error estimates using
Equation 16 of 0.15 years for going from 3 to 1 and 0.4 years for going from 1 to 3. In general, we found that the standard error was lower when the path follows the direction of flow. The top row of plots in Figure 8 appears to show that there is a slight bias, however the mean of the estimates are relatively close to our original estimate (original estimates are less than 0.3 years higher in both cases).5.2 Rotation
If we consider two points in the same H3 Index, for example location 1 () and a new point , then using the original grid system the method will simply produce 0 for this location as can be seen in Figure 9. We consider using 1000 rotations as explained in Section 4.1. For each rotation we estimate the travel time both back and forth. In 288 of the rotations the two points ended up in the same hexagon, hence resulting in a zero travel time. We plot the distribution of the other 712 travel times in the bottom row of Figure 9. The mean of all the entries including the zeros is 8.2 days for going from the new point to location 1, and 9 days for going from location 1 to the new point. This matches what we would expect as the currents in Figure 1 show no strong current going north or south near location , hence we would not expect a large difference in travel time in each direction.
The second benefit of performing rotations is to make the estimates less dependent on the grid system. We use the same 1000 rotations as with the previous example, and compute the most likely path and the mean travel times. In Figure 10 we plot the first 100 of the rotations along with the mean of the travel times resulting from these 100 rotations. The travel times and paths shown in this figure are comparable to those given in Figure 5. In most of the pathways we see that the darkest, most popular paths match up with the pathways in Figure 5.
One of the more interesting results from this analysis is the path going from to in Figure 10 a). Most of the paths go up closer to the Equator then use the Equatorial Counter current, followed by the Guinea and Gulf of Guinea currents as in the original vanilla application of the methodology. A small number of the rotations result in pathways that end up crossing the South Atlantic, to the south of location 2, then follows the South Atlantic current over to location 1. We expect this may be due to some of the grid cells blurring together two flow patterns, potentially a northward and southward current, resulting in the pathway along the West coast of Africa being less likely.
In general, the travel times from the rotation and original method can be quite different, which supports the need for this rotation methodology. If we compare Figure 5 and Figure 10, very few of the distances stay close to what they were in the original results using no rotations. In f) we see that going from 6 to 3 actually drops from 5.4 years to 3.1 years, with the most used path being very different to the original (no longer goes by location 2). This drop even causes the ordering of the distances to change as 6 to 3 is now the shorter travel time. Similarly, the ordering in e) changes, however the reasoning for this requires further investigation.
In Figure 6 we can see that the method does provide travel times from the South Atlantic into the Mediterranean Ocean despite there never being a drifter to cross the Strait of Gibraltar. In 64.6% of the 1000 rotations we were unable to obtain a travel time estimate from the Atlantic into the Mediterranean and in 96% we were unable to find a travel time estimate from the Mediterranean to the Atlantic. When we do not do a rotation we are able to obtain an estimate both in and out of the Mediterranean, this is due to the way the grid aligns as shown in Figure 3. Even if only one of the 1000 rotations are unable to provide an estimate it would be advisable to not use an estimate from this method. Therefore, using the vanilla method to estimate travel times into the Mediterranean is not a good option. Further adaptations to the method to provide added robustness to travel time estimates are discussed in Section 6.
6 Discussion and Conclusion
In contrast to van Sebille (2014), our methodology does not take into account seasonality. We have a few ideas for how seasonality could be incorporated. An obvious adaptation, if the aim was to obtain a short travel time which is expected to lie in a small 3 month window, is to just estimate using drifter observations which are in that window. Alternatively, we could use to be a certain jump such as a gap of two months, then we estimate 6 transition matrices say , where the entries are probabilities of going from the previous time period at state to state at the current time. Such a set up could still be solved using our shortest path algorithm. Other similar works such as Rypina et al. (2017) are more easily extendable to include seasonality. We justify our approach in the same way as Maximenko et al. (2012): we aim to provide a global view and a simple general concept explaining the pattern of potential pathways and travel times.
The use of the bootstrap and rotations are relatively easy methods to implement, each of which provides effective estimates of uncertainty. However, combining these procedures into one is not as simple. If we wanted to run rotations and bootstraps for each rotation, we still require a method to combine these estimates of travel times. We could treat every rotation equivalently, so say that our bootstrap sample in Equation 16 is all samples to obtain an estimate of uncertainty in travel time due to both grid discretization and data randomness. Additionally, we could decompose the uncertainty and provide a standard error for just the data randomness by estimating a standard error for each rotation using just the samples in each rotation, and then taking the average of all standard error estimates.
A better spatial indexing system may also yield further improvements. The H3 spatial index (UBER, 2019) is designed for land use, where the system is designed such that all the areas with highest distortion are placed within the ocean. The inequality in area is almost as bad as a longitudelatitude grid system, where the area of the largest hexagon is at most 1.6 times larger than area of the smallest hexagon. The adaptation to the grid described in Section 4.1 lessens the importance of this as the different sized grid cells will be placed in different areas under different rotations. Note that a regular longitudelatitude system could not be used with the rotations introduced in Section 5.2. The rotation would likely rotate the grids at a high latitude to an area of interest resulting in very small grid boxes, e.g. a grid at latitude could be around in area.
The method provided depends on the availability of drifter data making a connection at some point. Connections such as going across the Strait of Gibraltar are in practice impossible; any pathway which crosses it is due to a grid covering both the east and west of the Strait of Gibraltar. One potential way to adapt the method to approximate travel times across the Strait is, either adding artificial simulated trajectories as in van Sebille et al. (2012), or simply add a very small probability to the transition matrix crossing from the west to the east of the Strait of Gibraltar (and vice versa). For example, take two locations, one west and one east of the Strait of Gibraltar, say these correspond to states and respectively. If we wanted the crossing time to be 100 days into the Mediterranean set , the transition matrix will no longer be valid as the row no longer sums to one but the method will still work as intended. Such an adaptation would require the removal of the state which covers the Strait of Gibraltar to force the algorithm to take the artificial 100 day crossing.
The example with the Mediterranean Sea given in Section 5.2 is an interesting bonus feature of the rotation methodology but it is not as easily applicable to the Panama land mass problem. In the case of Panama, we will still obtain a travel time estimate from the Gulf of Mexico to the Pacific, but the times which are allowed to skip over the Panama land mass will be much shorter. An automatic detection could be achieved by looking at a large sample of rotations then running a test for multi modality. If it finds that there are two modes which are very far apart then this would be a sign that the method is finding some shortcut which is only present under some rotations. If such a method worked to detect the Panama land mass, we could then use it to search for more subtle surface transport barriers. In general it is preferable to preprocess the transition matrix such that rows/columns corresponding to unwanted links such at the Panama Canal and the Strait of Gibraltar are simply removed.
Our choice of the Lagrangian decorrelation time of 5 days may be too low in some instances. Previous works have found correlations in the velocity data lasting longer than days in certain regions (Lumpkin et al., 2002; Zhurbas and Oh, 2004; Elipot et al., 2010). This may suggest that using a larger value for may be needed to justify the Markov assumption. The tradeoff however is resolution, where shorter timescales allow pathways and distances to be computed with more detail. Our methodology is designed flexibly such that the practitioner can pick the most appropriate timescale for the spatial region and application of interest.
In general some unexpected features of the method do occur such as the discontinuity line in Figure 6 at approx latitude in the Pacific ocean. We expect there would be less of a discontinuity if these times were computed with the rotation methodology, however we argue that the discontinuities with travel times of most likely pathways should always exist. If smoothness of travel times was a major requirement, then one could consider the shortest path in travel time rather than the most likely path. The only necessary adaptation would be to use Equation 10 as the objective function in the shortest path algorithm rather than the negative log probability of Equation 7 that was used here. We expect the results would require more careful checking in such an approach, as the shortest path would be more likely to use any glitches in the grid system such as if there was a connection over Panama.
To summarize, in this paper we have created a novel method to estimate Lagrangian pathways and travel times between oceanic locations, thus offering a new, fast and intuitive tool to improve our knowledge of the dynamics of marine organisms and oceanic global circulation. Note that upon publication a Python package will be made available for public use which implements all the methodology presented.
References
 Arvo (1992) Arvo, J. (1992). Fast random rotation matrices. In KIRK, D., editor, Graphics Gems III, pages 117 – 120. Morgan Kaufmann, San Francisco.
 Becking et al. (2006) Becking, L. E., Cleary, D. F., de Voogd, N. J., Renema, W., de Beer, M., van Soest, R. W., and Hoeksema, B. W. (2006). Beta diversity of tropical marine benthic assemblages in the Spermonde Archipelago, Indonesia. Marine Ecology, 27(1):76–88.
 Berline et al. (2014) Berline, L. O., Rammou, A. M., Doglioli, A., Molcard, A., and Petrenko, A. (2014). A ConnectivityBased EcoRegionalization Method of the Mediterranean Sea. PLoS ONE, 9(11):1–9.
 Bonjean and Lagerloef (2002) Bonjean, F. and Lagerloef, G. S. E. (2002). Diagnostic model and analysis of the surface currents in the tropical pacific ocean. Journal of Physical Oceanography, 32(10):2938–2954.
 Dijkstra (1959) Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 1(1):269–271.
 Efron (1993) Efron, B. (1993). An introduction to the bootstrap. Monographs on statistics and applied probability ; 57. Chapman & Hall, New York.
 Elipot et al. (2016) Elipot, S., Lumpkin, R., Perez, R. C., Lilly, J. M., Early, J. J., and Sykulski, A. M. (2016). A global surface drifter data set at hourly resolution. Journal of Geophysical Research: Oceans, 121(5):2937–2966.
 Elipot et al. (2010) Elipot, S., Lumpkin, R., and Prieto, G. (2010). Modification of inertial oscillations by the mesoscale eddy field. Journal of Geophysical Research: Oceans, 115(9):1–20.
 Ellingsen and Gray (2002) Ellingsen, K. and Gray, J. (2002). Spatial patterns of benthic diversity: Is there a latitudinal gradient along the Norwegian continental shelf? Journal of Animal Ecology, 71:373 – 389.
 ESR (2009) ESR (2009). OSCAR third degree resolution ocean surface currents. Ver. 1. PO.DAAC, CA, USA. Dataset accessed [20200128] at https://doi.org/10.5067/OSCAR03D01.
 Gallo and Pallottino (1988) Gallo, G. and Pallottino, S. (1988). Shortest path algorithms. Annals of Operations Research, 13(1):1–79.
 Hansen and Poulain (1996) Hansen, D. V. and Poulain, P.M. (1996). Quality control and interpolations of wocetoga drifter data. Journal of Atmospheric and Oceanic Technology, 13(4):900–909.
 Huntley et al. (2011) Huntley, H. S., Lipphardt Jr, B., and Kirwan Jr, A. (2011). Lagrangian predictability assessed in the East China Sea. Ocean Modelling, 36(12):163–178.
 Jönsson and Watson (2016) Jönsson, B. F. and Watson, J. R. (2016). The timescales of global surfaceocean connectivity. Nature communications, 7(1):1–6.
 LaCasce (2008) LaCasce, J. H. (2008). Statistics from Lagrangian observations. Progress in Oceanography, 77(1):1–29.
 Lumpkin and Pazos (2007) Lumpkin, R. and Pazos, M. (2007). Measuring surface currents with surface velocity program drifters: the instrument, its data, and some recent results. Lagrangian analysis and prediction of coastal and ocean dynamics, 2:39–67.
 Lumpkin et al. (2002) Lumpkin, R., Treguier, A.M., and Speer, K. (2002). Lagrangian Eddy Scales in the Northern Atlantic Ocean. Journal of Physical Oceanography, 32(9):2425–2440.
 Maximenko et al. (2012) Maximenko, N., Hafner, J., and Niiler, P. (2012). Pathways of marine debris derived from trajectories of Lagrangian drifters. Marine Pollution Bulletin, 65(13):51–62.
 McAdam and van Sebille (2018) McAdam, R. and van Sebille, E. (2018). Surface connectivity and interocean exchanges from drifterbased transition matrices. Journal of Geophysical Research: Oceans, 123(1):514–532.
 Meehl (1982) Meehl, G. A. (1982). Characteristics of surface current flow inferred from a global ocean current data set. Journal of Physical Oceanography, 12(6):538–555.
 Miron et al. (2019) Miron, P., BeronVera, F. J., Olascoaga, M. J., and Koltai, P. (2019). Markovchaininspired search for MH370. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(4).
 Rypina et al. (2017) Rypina, I. I., Fertitta, D., Macdonald, A., Yoshida, S., and Jayne, S. (2017). Multiiteration approach to studying tracer spreading using drifter data. Journal of Physical Oceanography, 47(2):339–351.
 Smith et al. (2018) Smith, T. M., York, P. H., Broitman, B. R., Thiel, M., Hays, G. C., van Sebille, E., Putman, N. F., Macreadie, P. I., and Sherman, C. D. (2018). Rare longdistance dispersal of a marine angiosperm across the Pacific Ocean. Global Ecology and Biogeography, 27(4):487–496.
 Sykulski et al. (2016) Sykulski, A. M., Olhede, S. C., Lilly, J. M., and Danioux, E. (2016). Lagrangian time series models for ocean surface drifter trajectories. Journal of the Royal Statistical Society. Series C: Applied Statistics, 65(1):29–50.
 UBER (2019) UBER (2019). H3 spatial index. https://eng.uber.com/h3/. Accessed: 20200108.
 van Sebille (2014) van Sebille, E. (2014). Adrift.org.au  A free, quick and easy tool to quantitatively study planktonic surface drift in the global ocean. Journal of Experimental Marine Biology and Ecology, 461:317–322.
 van Sebille et al. (2011) van Sebille, E., Beal, L. M., and Johns, W. E. (2011). Advective time scales of Agulhas leakage to the North Atlantic in surface drifter observations and the 3D OFES model. Journal of Physical Oceanography, 41(5):1026–1034.
 van Sebille et al. (2012) van Sebille, E., England, M. H., and Froyland, G. (2012). Origin, dynamics and evolution of ocean garbage patches from observed surface drifters. Environmental Research Letters, 7(4).
 van Sebille et al. (2018) van Sebille, E., Griffies, S. M., Abernathey, R., Adams, T. P., Berloff, P., Biastoch, A., Blanke, B., Chassignet, E. P., Cheng, Y., Cotter, C. J., Deleersnijder, E., Döös, K., Drake, H. F., Drijfhout, S., Gary, S. F., Heemink, A. W., Kjellsson, J., Koszalka, I. M., Lange, M., Lique, C., MacGilchrist, G. A., Marsh, R., Adame, C. G. M., McAdam, R., Nencioli, F., Paris, C. B., Piggott, M. D., Polton, J. A., Rühs, S., Shah, S. H. A. M., Thomas, M. D., Wang, J., Wolfram, P. J., Zanna, L., and Zika, J. D. (2018). Lagrangian ocean analysis: Fundamentals and practices. Ocean Modelling, 121:49–75.
 van Sebille et al. (2009) van Sebille, E., van Leeuwen, P., Biastoch, A., Barron, C., and de Ruijter, W. (2009). Lagrangian validation of numerical drifter trajectories using drifting buoys: Application to the Agulhas system. Ocean Modelling, 29(4):269–276.
 Wakata and Sugimori (1990) Wakata, Y. and Sugimori, Y. (1990). Lagrangian motions and global density distributions of floating matter in the ocean simulated using shipdrift data. Journal of Physical oceanography, 20:125–138.
 White et al. (2010) White, C., Selkoe, K. A., Watson, J., Siegel, D. A., Zacherl, D. C., and Toonen, R. J. (2010). Ocean currents help explain population genetic structure. Proceedings of the Royal Society B: Biological Sciences, 277(1688):1685–1694.
 Zhurbas and Oh (2004) Zhurbas, V. and Oh, I. S. (2004). Drifterderived maps of lateral diffusivity in the Pacific and Atlantic oceans in relation to surface circulation patterns. Journal of Geophysical Research: Oceans, 109(C5).
Comments
There are no comments yet.