We used transition path theory (TPT) to infer “reactive” pathways of floating marine debris trajectories. The TPT analysis was applied on a pollution-aware time-homogeneous Markov chain model constructed from trajectories produced by satellite-tracked undrogued buoys from the National Oceanic and Atmospheric Administration's Global Drifter Program. The latter involved coping with the openness of the system in physical space, which further required an adaptation of the standard TPT setting. Directly connecting pollution sources along coastlines with garbage patches of varied strengths, the unveiled reactive pollution routes represent alternative targets for ocean cleanup efforts. Among our specific findings we highlight: constraining a highly probable pollution source for the Great Pacific garbage patch; characterizing the weakness of the Indian Ocean gyre as a trap for plastic waste; and unveiling a tendency of the subtropical gyres to export garbage toward the coastlines rather than to other gyres in the event of anomalously intense winds.

Given a Markov chain, namely, a model describing the stochastic state transitions in which the transition probability of each state depends only on the state attained in the previous event, transition path theory (TPT) provides a rigorous approach to study the statistics of transitions from a set of states to another, possibly disconnected set of states. Envisioning the motion of floating debris as described by a Markov chain that accounts for the ability of coastal states to “pollute the oceans,” TPT is employed to unveil “reactive” pathways representing direct transitions from potential release locations along the shorelines to accumulation sites across the world ocean. These include the subtropical gyres, whose strength in this context is investigated.

The long-term fate of satellite-tracked drifting buoys from the National Oceanic and Atmospheric Administration's (NOAA) Global Drifter Program1 is characterized by a tendency to form clusters in the oceans’ subtropical gyres2,3 that resemble great garbage patches.4 The development of such clusters, most evidently in the case of undrogued (i.e., without a sea anchor) drifters,5 has been explained5–7 as the result of the combined action on the drifters of converging ocean currents and winds mediated by their inertia, which prevent them from adapting their velocities to that of the carrying water–air flow system.

The tendency of the drifters to cluster in the long run enables a probabilistic description of their dynamics using results from ergodic theory8 and Markov chains,9,10 which form the basis for approximating asymptotically invariant sets using so-called set-oriented methods11–14 and related methods including flow networks.15–17 This approach places the focus on the evolution of probability densities, which, unlike individual trajectories, represent robust features of the dynamics. Central to this measure-theoretic characterization is the transfer operator and the transition matrix, its discrete version resulting by covering the phase space with boxes, which represent the states of the associated Markov chain.

Such a probabilistic description has been applied on simulated drifter trajectories,18 suggesting a characterization of great garbage patches as almost-invariant attracting sets with corresponding basins of attraction spanning areas as large as those of the geographic ocean basins. While the latter suggests a strong influence of the regions collecting marine debris on their global transport, it does not provide information on pollution routes.

The goal of this paper is to unveil such routes from observed drifter trajectories. This is done by applying transition path theory (TPT).19–23 Developed to investigate transition pathways in complex nonlinear stochastic systems, TPT provides a statistical characterization of the ensemble of “reactive” trajectories, namely, pieces of trajectories along which direct transitions between two sets A and B in phase space take place. The TPT terminology is borrowed from statistical mechanics and physical chemistry, for which TPT was originally developed to study chemical reactions from reactants A to products B, as an improvement for earlier approaches such as transition state theory24 and transition path sampling.25 Since then, the TPT framework has also been applied to study molecular conformation changes26,27 and transitions in climate models.28,29 We here present, to the best of our knowledge, the first oceanographic application.

By constructing a Markov chain for debris motion and then identifying coastline boxes in the ocean covering with reactant states A, and boxes in several ocean locations including the subtropical gyres with product states B, we use TPT to infer pollution pathways in the global ocean. The Markov chain model accounts for the ability of coastal boxes (states) to “pollute the oceans.” This involves adding an artificial state to the chain where all outflow goes in and all inflow comes from (in a manner that differs from prior approaches30,31). By setting A to a single garbage patch and B as the union of the other garbage patches, we can also assess the strength of the patches.

The rest of the paper is organized as follows. The ergodic-theory setup for closed systems is presented in Sec. II. An adaptation of the theory for open systems in discussed in Sec. III. The main results of TPT are reviewed in Sec. IV, both for closed systems (Sec. IV A) and an extension for open domains (Sec. IV B). The Markov-chain model for ocean pollution is constructed in Sec. V from satellite-tracked drifter trajectories. This entails coping with a number of issues, previously not encountered, partially addressed, or overlooked;32–36 these include zonal connectivity, spurious communication between ocean basins, and nonobserved communication, as well as incorporating pollution sources near the coast. In Sec. VI, time-asymptotic aspects of the chain dynamics are investigated, suggesting prospects for garbage patches yet to be directly observed. The TPT analysis is applied in Sec. VII. This reveals pollution routes into the garbage patches, which represent alternative targets for ocean cleanup efforts.37 Finally, a summary and the conclusions of the paper are presented in Sec. VIII.

Let us assume that floating debris trajectories are described by a time-homogeneous stochastic process in continuous space XR2 and observed at discrete times nT, nZ. Its transition probabilities are controlled by a stochastic kernel K(x,y)0 such that XK(x,y)dy=1 for all x in phase space X, representing the world ocean basin. The stochastic kernel is time-independent since the time-homogeneity of the process implies that the rules governing the process at any time are the same. It is convenient to think of X as a measure space, i.e., a set equipped with a σ-algebra of subsets measured by (normalized) area. Then, a probability density f(x)0, Xf(x)dx=1, describing the distribution of the random position XnT at any time nT evolves to the distribution

(1)

at time (n+1)T, which defines a Markov operator P:L1(X) generally known as a transfer operator.8 

To infer the action of P on a discretized space, one can use a Galerkin projection referred to as Ulam’s method.14,38,39 This consists of covering the phase space X with N connected boxes {Bi}iS, S:={1,,N}Z+, disjoint up to zero-measure intersections, and projecting functions in L1(X) onto the finite-dimensional space spanned by indicator functions on the boxes VN:=span{1Bi(x)area(Bi)}iS, where 1A(x)=1 if xA and 0 otherwise. The discrete action of P on VN is described by a matrix P=(Pij)i,jSRN×N called a transition matrix. The transition matrix results from the projection,33,34

(2)

describes the proportion of probability mass in Bi that flows to Bj during T. If one is provided with a large set of observations x0 and xT of X0 and XT, respectively, then (2) can be estimated via counting the transitions in the observed data, viz.,

(3)

Note that jSPij=1 for all iS, so P is a row-stochastic matrix that defines a Markov chain on boxes, which represent the states of the chain.9,10 The evolution of the discrete representation of f(x), i.e., the probability vector f=(fi)iS, iSfi=1, is calculated under left multiplication, i.e.,

(4)

as it follows by noting that Pr(X(n+1)TBj)=iSPr(X(n+1)TBj,XnTBi)=iSPr(XnTBi)Pij. In this paper, whenever we multiply vectors by matrices, we assume that the vector takes the appropriate form of a row or column vector for the given operation.

Because P is stochastic, 1=(1,,1) is a right eigenvector with eigenvalue λ=1, i.e., P1=1. The eigenvalue λ=1 is the largest eigenvalue of P. The associated potentially nonunique left eigenvector p=(pi)iS is invariant, because pP=p and can be chosen componentwise nonnegative (by the Perron–Frobenius theorem).

We call Pirreducible (or ergodic) if for all i,jS there exists nijZ0+{} such that (Pnij)ij>0. To wit, all states of an irreducible Markov chain communicate, the eigenvalue λ=1 is simple, and the corresponding left eigenvector p is strictly positive.10 We call Paperiodic (or mixing) if there exists iS such that gcd{nZ0+:(Pn)ii>0}=1. No state of an aperiodic Markov chain is visited cyclically.

If P is ergodic and mixing, then p, normalized to a probability vector (iSpi=1), satisfies 0<p=pP=limnfPn for any probability vector f. We call p an invariant limiting probability vector or stationary distribution.

We adopt the traditional notation with {Xt}tZ instead of {XnT}nZ and write, for instance, Pij=Pr(Xt+1=jXt=i), when this simplifies the notation. In what follows, we will assume that P is both ergodic and mixing, and the system is in stationarity, i.e., Pr(XtBi)=pi for all tZ.

The Markov chain model we will deduce from data in Sec. V A is, however, open, thus not ergodic. For this reason, we shall next consider the closure of open dynamics.

Let us assume that the flow domain is no longer closed, meaning that trajectories can flow out of the domain and back into it. This can happen, for instance, when the domain of interest is a subregion of the closed world ocean domain X or when trajectory data are only available in a subregion of X. Other possibilities include poor sampling of X, weak communication within, or the situation we describe in Sec. V. In every case, the resulting dynamical system represents an open dynamical system.

The above is a slight variation of the setting in Sec. II. We still assume that the motion is described by a discrete-time-homogeneous Markov chain on a box covering {Bi}iO of the ocean domain X but the probability to transition from one box with index iO to anywhere else in the domain O is no longer strictly 1 since probability mass can flow out of the domain. We denote the transition matrix on the open domain by PO with entries given by PijO:=Pr(Xt+1=jXt=i) for i,jO. Since the rows of PO no longer have to add up to one, PO represents a substochastic matrix.

We assume that a larger domain SO exists on which the dynamics are closed, i.e., the transition matrix P on box entries i,jS is stochastic. Furthermore, when we say that the dynamics on the open domain is stationary, we actually mean that the dynamics on the larger, closed domain is stationary with distribution p=(pi)iS, while we denote the restriction to the open domain by p|O=(pi)iO.

For further analysis, it is often useful to artificially close the open system. From the closure of PO, we can, for instance, get an estimate of p|O. Closing PO can be done by appending to O a state ω, which we will call a two-way nirvana state, and letting all the outflow from O flow into ω, while also redistributing the probability mass from ωback into O. Since thereby all boxes that are in S but not in O are lumped together, this restricted dynamics should be consistent with the original one under the assumption of well-mixedness between exit from O and reentry into it. For simplicity of notation, we will denote the singleton {ω} also by ω and refer to it too as the two-way nirvana state.

The resulting transition matrix on Oω reads (possibly overloading the notation by denoting it by P again)

(5)

where POω:=(1jOPijO)iO (understood as a column vector) gives the outflow from O to ω and PωO is a (row) vector that gives the inflow and has to be a probability vector. Note that the matrix P is stochastic jOωPij=1 for all iOω and as such constitutes a closed dynamical system.

When no information about the reentry is available, e.g., because data outside the open domain of interest are not available, a possible choice30 for PωO is to redistribute according to the quasistationary distribution of PO. Lünsmann and Kantz31 alternatively use contour advection for estimating the transition probabilities between boxes. Without adding a nirvana state, Froyland et al.30 immediately redistribute the outflow back into the system. Here, we redistribute in such a way that accounts for ocean pollution, as we describe in Sec. V.

In Sec. IV, we will see how to study transitions between A and B (subsets of O) in both the cases where (i) the domain is closed, i.e., O=S, and where (ii) paths only traverse the open domain OS. In the latter case, for the TPT computations, only knowledge of PO and the estimate of the stationary density on the open computational domain p|O is necessary.

Motivated by a desire to understand rare events such as transformations involved in chemical reactions, TPT provides a rigorous approach to study transitions from a set AS to another, disjoint set BS of a Markov chain. The results presented below pertain to time-homogeneous (i.e., autonomous) chains;19–22 extensions to the nonautonomous case have been recently derived,40 but they are beyond the scope of this paper. Traditionally, source set A is thought to be formed by reactant states, while target set B of product states. Thus, transitions from A to B are referred to as reaction events, while the pieces of trajectories running from A to B without going back to A or going through B in between are known as reactive trajectories, which are the focus of TPT (Fig. 1).

FIG. 1.

Given a Markov chain taking values on S, the cartoon shows in red the reactive pieces of a trajectory connecting disjoint sets A,BS.

FIG. 1.

Given a Markov chain taking values on S, the cartoon shows in red the reactive pieces of a trajectory connecting disjoint sets A,BS.

Close modal

The main tools of TPT are the forward and backward committor probabilities giving the probability of a random walker to hit B before A, in either forward or backward time. The committor probabilities are used to express various statistics of the ensemble of reactive trajectories: (i) the density of reactive trajectories, which provides information about the bottlenecks during the transitions; (ii) the current of reactive trajectories indicating the most likely transition channels; (iii) the rate of reactive trajectories leaving A or entering B; and (iv) the mean duration of reactive trajectories. We will introduce these in the following. Recall that we assume the chain to be stationary with distribution p.

The first entrance time of a set SS is the stopping time random variable defined as

(6)

where inf:=. The forward committorq+:=(qi+)iS gives the probability that a trajectory starting in iSfirst entersB, notA, i.e.,

(7)

Note that qiA+=0 while qiB+=1. For iC:=S(AB), one has that

(8)

The solution to this algebraic system is unique due to the irreducibility of P, and in matrix notation expressed as

(9)

where |S denotes the restriction on indices in S, while |S,S gives the restriction to rows corresponding to S and columns of S; if S=S, we shorten this to |S.

The last exit time, in turn, is defined by

(10)

where sup:=, which is a stopping time, but for the time-reversed chain {Xt}tZ that traverses the original Markov chain backward in time, i.e., Xt:=Xt. The reversed chain’s transition matrix, P=(Pij)i,jS, is given by

(11)

since the chain is assumed to be in stationarity. The time-reversed transition matrix P is ergodic and mixing and has the same stationary distribution p as P. The backward committorq:=(qi)iS gives the probability that a trajectory starting in iSlast exitsA, notB,

(12)

In this case,

(13)

for iC, subject to qiB=0 and qiA=1. The (unique) solution in matrix notation,

(14)

A particular situation arises in the special case when the chain is reversible, namely, when piPij=pjPji or, equivalently, P=P. In such a case, q=1q+.

The committors contain information that enable the computation of various transition statistics. The distribution of reactive trajectoriesμAB=(μiAB)iS, defined as the joint probability that the chain is in state i while transitioning fromAtoB, viz.,

(15)

tells us where reactive trajectories spend most of their time. Note that μiABAB=0. The distribution of reactive trajectories is computable from the committor probabilities and the stationary distribution,

(16)

A density of reactive trajectoriesμ^AB=(μ^iAB)iS is obtained by normalizing μiAB by the probability to be reactive,

(17)

as it follows from the law of total probability. The result is

(18)

i.e., the probability of being in state i conditioned on being already on a reactive path from A to B.

The current (or flux) of reactive trajectoriesfAB=(fijAB)i,jS gives the average flux of trajectories going through i and j at two consecutive times while on their way from A to B,

(19)

which is computable as

(20)

Note that the reactive current can include direct transitions from iA to jB, which are not accounted for in the corresponding reactive distribution as it only considers transitions passing through C.

To eliminate detours of reactive currents, one introduces the effective current of reactive trajectoriesf+=(fij+)i,jS, which gives the net amount of reactive current going through i and j consecutively, viz.,

(21)

To visualize f+ on a flow domain covered by boxes {Bi}iS, one usually depicts the magnitude and the direction of the effective current out of each i, i.e., to each i one attaches the vector jifij+eij, where eij is the unit vector pointing from the center of box Bi to the center of Bj. There also exists a flow decomposition algorithm for extracting the dominant transition paths from f+.22 

The rate of transitions leaving A or departure rate is defined as the probability per time step of a reactive trajectory to leave A, i.e.,

(22)

and can be computed by summing up the reactive flux that exits A. In turn, the rate of transitions entering B or arrival rate is defined as the probability per time step of a reactive trajectory to enter B,

(23)

By a simple calculation, it can be shown that summing the reactive current out of A, iA,jSfijAB, is equal to aggregating the reactive current into B, iS,jBfijAB, thus

(24)

To better interpret the transition rate kAB, we give two meanings. Consider an infinite p-distributed ensemble of random walkers in our domain, then at any time the proportion of random walkers that are exiting A while on their way to B (or equivalently, entering B when coming last from A) is given by kAB. Now, on the other hand, consider only one random walker in the system, then kAB can be interpreted as a frequency, i.e., the random walker exits A on average every (kAB)1-th time on the way to B (and, equivalently, enters B when coming from A).

In some situations, e.g., when B is given by a disconnected set, it is insightful to further decompose the transition rate

(25)

into the individual arrival rates into disjoint subsets Bn that together give B=nBn,

(26)

The same can also be done for decomposing kA.

Finally, dividing the probability of being reactive by the discrete transition rate,

(27)

gives the expected duration of a transition from A to B.19,40

We close this section with a remark on comparing probabilistic computations with counting. Ergodicity of the chain implies that the objects in TPT can be approximated by “counting” transition events of one sufficiently long trajectory, and this approximation converges almost surely as the length of the trajectory tends to infinity.19,40 For instance, the forward committor qi+ of any state i is approximated by the fraction of all visits of the chain to state i after which the chain directly transitioned to B without hitting A first. All other quantities considered here can be similarly approximated. As we intend to apply TPT to a chain extracted from drifter trajectory data, one might wonder whether this level of sophistication is necessary to our ends or whether one could simply do an approximation by counting. The answer lies in the features of the data. One would need sufficiently many drifter trajectories that are sufficiently long to resolve the transition statistics and that are also spread according to the right distribution. None of these requirements are met, and the best one can do is to “concatenate” the drifter information into a Markov chain, as it will be done in Sec. V below.

To apply TPT to open dynamical systems on O, a modification from the standard setting as reviewed in Sec. IV A is needed. Adding the state ω to O closes the system artificially (as in Sec. III), but we are still only interested in the transitions from AO to BO that stay in O during the transition. Thus, the reactive trajectories we consider go from A to B without passing A, B, or ω during transition. If we were to apply the usual TPT on the artificially closed system, we would also observe artificial transitions via the added state ω.

In order to compute the statistics of the reactive trajectories from A to B only through O, we look at slightly different committors. Namely, the forward committor now gives the probability to next transition to B rather than to A or outside of O when starting in state i, i.e.,

(28)

while the backward committor gives the probability to have last come from A, not Bω,

(29)

In that way, the product of forward and backward committors becomes the probability when initially in i to have last come from A and next go to B while not passing through A, B, or ω in between.

By definition, the forward committor is qi+=0 for iAω and 1 for iB, while in the transition region C:=O(AB), it satisfies

(30)

since qω+=0 and P on entries of O reduces to PO.

The backward committor qi=0 for iB and 1 for iAω, while, by a similar reasoning as above, it satisfies

(31)

for iC, where PO, is the restriction of the backward-in-time transition matrix P to O and has entries PijO,=pjpiPjiO for i,jO.

Therefore, system (9) remains the same with the replacement of P with PO and A with Aω. In turn, system (14) remains the same with the replacement of P with PO, and B with Bω.

The rest of the formulae in Sec. IV A are not changed except that the committors are now given as above. An important observation, however, is that μiAB=0 for i=ω and fijAB=0 for i,j=ω. Thus, only their values on O are of interest, where P can be replaced by PO and p can be substituted by its restriction to O, p|O. Also, as the rate and mean transition time of reactive trajectories are derived from density and current, they are computable solely from PO and p|O.

This version of TPT for open dynamics can also be useful in other settings, e.g., when one wants to study transitions between A and B that avoid a third subset D of the state space S.

In the following, we describe our stochastic model for the dynamics of a single plastic debris piece that enters the ocean at the coast with a probability reflecting observed levels of mismanaged plastic waste in near coastal communities. From the coast, the debris piece traverses the ocean, possibly passing and staying for long times near garbage patches. Its motion is fitted using satellite-tracked drifter trajectories; cf. Sec. V A. Whenever a debris piece beaches somewhere, we reinject it again next to the coast. The coastal injection and beaching are described in Sec. V B below. Note that the results of TPT analysis depend on the underlying distribution of the system. Without reinjection, we would be dealing with an open system, and it is not obvious what distribution would be natural to consider. We here choose the saturated plastic distribution that results from the constant injection rate and beaching probability in the infinite time limit. It turns out that reinjecting the beached process closes the otherwise open system in such a way that the stationary closed system is equivalent to the saturated open one. This is shown in Sec. V C. Any other choice would require a non-stationary TPT analysis, which is beyond the scope of this article.

As anticipated, to formulate the Markov chain for marine debris motion, we use drifter trajectory data, taken from the NOAA Global Drifter Program.1 Satellite-tracked by the Argos system or GPS (Global Positioning System), the drifters from this database have a spherical surface float with a 15-m-long holey-sock drogue attached.41 They are engineered to resist wind slippage and wave-induced drift and hence to follow water motion as close as possible.42 We therefore only consider trajectory portions during which the drifter’s drogue has been lost,43 which can be expected to provide a more fair representation of floating marine debris motion.5–7,44,45

The basic procedure to construct the transition matrix P, defined in (3), is as follows. We first interpolate the available undrogued drifter trajectories daily and form two arrays, one representing positions at any instant of time over 1992–2019 (x0) and another one representing their images (xT) after T=5 d. Here, we are assuming that the ocean motion did not change considerably over the last 30 years such that the transition matrix P from this dataset is still a good representation of the “average ocean motion.” In other words, effects associated with climatic variability are ignored.

We then define the box covering {Bi}iS by lying down on the world ocean domain a grid of roughly 3° width (due the planet’s curvature the area of the boxes is not fixed, varying from 100 to 10 000 km2, but this is inconsequential in the definition of the vector space VN, normalized by box area). The entries of P are finally estimated via counting according to (3). As in previous work,32–36 the transition time T is chosen long enough to guarantee negligible memory into the past and sufficient communication among boxes, made large enough to maximize sampling. Even if the full system dynamics is described by a memoryless advection-diffusion process, when projected on a finite-state chain, it cannot in general be expected to result in a Markov process due to possible partition box boundary recrossings.46 In our case, the simple Markovianity test λ(P(nT))λ(P(T))n is passed well up to n=4, where λ() denotes the spectrum of the argument and P(t) is the estimated transition matrix for transition time t. This gives us confidence in our analysis.

There are additional aspects, not encountered, partially addressed, or overlooked earlier, which must be coped with to make P meaningful.

  1. Zonal connectivity. This is addressed by identifying and continuating trajectories crossing the antimeridian connecting the eastern and western hemispheres.

  2. Spurious communication between ocean basins. This situation occurs where ocean basins are separated by narrow land masses. The situations that concern us are the Panama Isthmus separating the Pacific and Atlantic Ocean basins, and also the Maritime Continent separating the Pacific and Indian Oceans. Neither the undrogued drifters considered nor drogued drifters analyzed earlier47 reveal connectivity between the Pacific and Indian Oceans through the various straits and passages in that region, which might seem at odds with the presence of the Indonesian throughflow,48 particularly for the drifters drogued at 15 m. However, this takes place mainly within the thermocline layer (50–200 m),49 which is less correlated with the local wind flow that quite strongly affects the undrogued drifters and also the drogued drifters, albeit to a lesser extent. To avoid spurious communication between the basins, we proceed as follows. Let Bk be a box spanning portions of, for instance, the Pacific Ocean and Atlantic Ocean (Caribbean Sea). Denote BkPO and BkAO the portions of Bk lying on the Pacific Ocean and Atlantic Ocean sides, respectively. In computing transitions between Bk and other boxes, we only consider those from or into BkPO or BkAO depending on which one makes the largest number of transitions. This guarantees that Pkj>0 and Pjk>0 exclusively for jS such that Bj is either in the Pacific Ocean or in the Atlantic Ocean.

  3. Nonobserved communication. A prominent example of this is the communication between the Atlantic Ocean and the Mediterranean Sea. Depending on the size of the boxes Bi, a connection might exists through the Gibraltar Strait, even thought in reality no drifter is seen to traverse it (in any direction). We resolve this situation by excluding the Mediterranean Sea domain from consideration.

  4. Weak communication. We enable as much communication as possible along the chain by restricting the chain to the largest strongly communicating class of states. This is done by applying the Tarjan algorithm50 on the directed graph equivalent to the Markov chain. This procedure excludes boxes from the partition. Among those boxes are 22 poorly sampled coastal boxes, mainly in the Kara Sea of the Arctic Ocean and the Seas of Indonesia, with trajectories flowing in, but not flowing out in the next step. Let O be the ordered set of box indices in the largest class of strongly connected boxes. Using the notation in Sec. III, we call PO the substochastic transition matrix characterizing this open system. The Markov chain is now substochastic, since by the exclusion of boxes, it is no longer ensured that probability mass is conserved.

To formulate our Markov-chain model for ocean pollution, we leverage the possibility that marine debris get stuck at shorelines. This creates an additional outflow of the system that must be compensated for, which we choose to do in such a way as to model ocean pollution at the coasts.

Specifically, let :O[0,1) be a land fraction function giving the ratio between the land area and the total box area. Namely, 0<(i)<1 for iLO corresponding to boxes filled with some portion of land (or ice) (Fig. 2, top panel) and (i)=0 otherwise. We then follow Miron et al.33 and replace PijO with

(32)

for all i,jO. To wit, only a fraction of the probability mass, proportional to the amount of land covering box Bi, is allowed to flow from i to j, the remaining probability mass is assumed to beach and flows out of the system. The factor α, not considered in Miron et al.,33 was included to enable consistency with observations. While we have performed optimizations of no kind, we have found that α=14 produces results most consistent with them. If α=1 (as in Miron et al.,33) then the so-called Great Pacific garbage patch51 in the North Pacific subtropical gyre is not revealed as intense as observations indicate.4,52 However, transition channels into this patch and patches in the other subtropical gyres are not sensitive to the specific α-value assumed, as we show in the supplementary material.

FIG. 2.

(Top panel) Fraction of land (or ice) filling coastal boxes of the surface world ocean partition (black). (Bottom panel) Percentage of share of global mismanaged littered or inadequately disposed plastic waste estimated in 2010 for populations living within 50 km of the coastline.

FIG. 2.

(Top panel) Fraction of land (or ice) filling coastal boxes of the surface world ocean partition (black). (Bottom panel) Percentage of share of global mismanaged littered or inadequately disposed plastic waste estimated in 2010 for populations living within 50 km of the coastline.

Close modal

To deal with the created substochasticity by a closure of the system, we augment the chain by one artificial state ω as in (5). All the outflow of the open system goes into ω and we reinject the probability mass from ω to O through coastal boxes according to the plastic waste input from land into the ocean, viz.,

(33)

Here, Wi is the mass of mismanaged plastic waste in Bi, iL, as inferred from estimates53 made in 2010 for populations living within 50 km of the coastline. This is shown in percentage of the total mass in the bottom panel of Fig. 2; note that only inhabited coastal boxes for which estimates are available are shown. We denote the transition matrix on the closed domain by P, but it should not be confused with the transition matrix P from Sec. V A, which has a different domain and entries. The two-way nirvana state ω compensates for the substochasticity of PO=P|O by sending back into the chain any imbalances through the land states distributed according to the ability of such states to “pollute the oceans” as inferred by their share of the global land-based plastic waste entering the ocean through them. It must be realized that in this statistical model debris mass is neither created nor destroyed. In other words, the model assumes that the world ocean is polluted by plastic at a certain level, and that the ocean currents and winds redistribute the existing pollutants within ocean basins. If beaching occurs, then the pollutants are returned back to the ocean in an equal quantity simulating mismanaged plastic waste loading from land runoff.

We model the distribution of garbage input per time unit to the oceans by a time-independent vector PωO=:rR|O|. Each entry of r accounts for the probability per time unit of injecting a garbage particle into the corresponding box. Thus, r is supported on the coastal (land) boxes, i.e., ri=0 for iL, and r is a probability vector, i.e., iri=1.

Then, the total accumulated garbage mass distribution in the oceans is time-asymptotically going to be

(34)

Recall that PO is assumed to be irreducible; thus, IdPO is invertible. Equation (34) gives the mass distribution of debris particles entered over an infinite time frame; thus, it does not need to be a probability vector. It is the limiting (saturated) mass distribution of pollution measured in the units dictated by r. If we would like to know the relative distribution of garbage that has accumulated over time, we would normalize this vector to a probability vector.

Now, it turns out that the very same long-term distribution is modeled by our “recirculating” Markov chain. With a:=POω being the vector of absorption probabilities from the boxes into nirvana (the outflow from the open domain), the stationary distribution of our chain satisfies

(35)

with stationary vector p|O on indices of O and scalar ρ giving the stationary weight of ω. The set of equations corresponding to the boxes in O read as p|OPO+ρr=p|O, or, after rearrangement,

(36)

Since ρ is scalar, this readily means that p|Or(IdPO)1, which equals the asymptotic mass distribution (34) from above. In summary, our stationary Markov chain constructed with reinjection is the statistical equivalent of garbage motion in the ocean based on the limiting garbage distribution.

By design, the proposed transition matrix P for marine debris pollution has a single maximal communicating class of the states, implying irreducibility for P and ergodicity for the dynamics. Furthermore, direct pushforward (i.e., evolution under left multiplication by P) of an arbitrary probability vector reveals convergence to the dominant left eigenvector p (the chain is also aperiodic), which is invariant and also limiting, and hence represents a stationary distribution. The top panel of Fig. 3 shows p>0 restricted to O, viz., the set of boxes of the world ocean partition where the dynamics are open. The middle and bottom panels show, restricted to O, the distribution after 1 and 10 years of evolution under left multiplication by P of 1ω, respectively. Note that p|O locally maximizes in the subtropical gyres, quite evidently in the eastern side of the North Pacific gyre. In most of the Indian Ocean, p|O reveals several well-spread local maxima consistent with a predominantly uniform distribution. The exception is the Bay of Bengal, where p|O shows more clear sings of local maximization. An additional local maximum of p|O is seen in the Gulf of Guinea south of West Africa. The several local maxima of p|O identified are indicated by the red boxes in the top panel of Fig. 3. The Indian Ocean location corresponds to its local maximum inside the subtropical gyre.

FIG. 3.

(Top panel) Restricted to the set O of boxes covering the physical ocean domain where the dynamics are open, the stationary distribution p of the closed dynamics represented by the transition matrix P for marine debris pollution. Note that p locally maximizes inside the subtropical gyres, which, at the same time happen to develop great patches. Indicated by the red boxes are these (and additional; cf. text for details) local maxima of p|O. (Middle panel) Restricted to O, distribution after 1-year evolution under left multiplication by P of a probability density (vector) with support on the virtual nirvana state included to close the system. (Bottom panel) As in the middle panel, but after 10 years.

FIG. 3.

(Top panel) Restricted to the set O of boxes covering the physical ocean domain where the dynamics are open, the stationary distribution p of the closed dynamics represented by the transition matrix P for marine debris pollution. Note that p locally maximizes inside the subtropical gyres, which, at the same time happen to develop great patches. Indicated by the red boxes are these (and additional; cf. text for details) local maxima of p|O. (Middle panel) Restricted to O, distribution after 1-year evolution under left multiplication by P of a probability density (vector) with support on the virtual nirvana state included to close the system. (Bottom panel) As in the middle panel, but after 10 years.

Close modal

The structure of p|O suggests garbage patches in the subtropical gyres of the Atlantic and Pacific Oceans consistent with in situ microplastic concentration observations.4 Previous analyses2,3 of drifter data revealed these patches too, albeit from direct evolution of probability densities. In particular, van Sebille et al.2 argued that the North Pacific patch should be the main attractor of global marine debris, in agreement with direct observational evidence4,52 of the Great Pacific garbage patch.51 Our pollution-aware model produces consistent results. This can be anticipated from p|O acquiring larger values in the North Pacific gyre than in the other subtropical gyres, and also from direct pushforward of 1ω and subsequent restriction of the evolved density to the boxes where p|O locally maximizes in the subtropical gyres (Fig. 4). (It should be noted too that the structure of p|O in the North Pacific suggests a garbage patch, albeit weaker, in the western side of the basin in agreement with field sampling.54) The relative weakness of the Indian Ocean garbage patch4 attributed to unique oceanic and atmospheric dynamics in the region55 is consistent with the results from our Markov-chain model for ocean pollution too. There the stationary distribution p|O does not reveal a clear local maximum (Fig. 3), and the direct pushforward of 1ω identifies the Indian Ocean gyre as the less attracting of all the subtropical gyres (Fig. 4). Exactly where the garbage patches are located is determined by wind-induced Ekman and wave-induced Stokes drift effects3 mediated by the inertia (i.e., buoyancy and size) of the floating debris pieces.5–7,44,45,56 Indeed, the numerical simulations of inertial particles by Beron-Vera et al.5 do not reveal signs of accumulation in the Indian Ocean gyre as clear as in the other gyres.

FIG. 4.

Evolution of 1ω under left multiplication by P restricted to the boxes where p|O locally maximizes in the subtropical gyres.

FIG. 4.

Evolution of 1ω under left multiplication by P restricted to the boxes where p|O locally maximizes in the subtropical gyres.

Close modal

In Sebille et al.,2 the authors suggest the possibility of a rather strong garbage patch in the Barents Sea in the Arctic Ocean, possibly constrained by slow surface convergence due to deep-water formation. While the authors noted that this patch might be an artifact of drifters becoming grounded in the (seasonal) sea-ice, observational support of plastic accumulating in the region is emerging.57 However, the observed accumulation represents a very small fraction (3%) of the global standing stock. Our pollution-aware model does not reveal a patch there, more consistent with this observation.

Our model suggests the occurrence of a patch in the Gulf of Guinea, which seems to be supported only on numerical simulations.58 However, the Gulf of Guinea is identified as a mesopelagic niche with genomic characteristics than different than its surroundings.59 This patch remained elusive to earlier studies.2,3 A likely explanation is the involvement in those earlier studies of both undrogued and drogued drifters, which unlike floating debris, are much less affected by inertial effects.56 However, a more recent study55 involving exclusively undrogued drifters did not reveal accumulation in the Gulf of Guinea time asymptotically.

The structure of p|O also reveals that the Bay of Bengal has potential for holding a garbage patch. High plastic concentration in the Bay of Bengal has been reported and attributed to loading from nearby land-based sources.60 The occurrence of a garbage patch in the Bay of Bengal was also suggested recently from the analysis of undrogued drifter trajectory data.55 

The pertinent question is how the garbage patches are filled. We address this using TPT.

With the above in mind, we proceed to apply TPT to the dynamics on the physical world ocean domain, where reactive debris currents are sought to be unveiled. The usual TPT (Sec. IV A) allows us to compute statistics of the ensemble of reactive paths of marine debris into garbage patches. With the help of TPT for open domains (Sec. IV B), we can study reactive paths between garbage patches.

To infer the pollution paths into garbage patches, we choose the nirvana state ω as the source state A of garbage, and we identify the set of target states B with the union of indices of boxes covering garbage patches as inferred by the regions where p|O tends to locally maximize, which we have isolated above (cf. Fig. 3, top panel). We will denote G the garbage patch set. Although the debris is reentering the ocean through the land boxes L, choosing the source as A=ω is more reasonable, as it allows reactive debris trajectories to enter boxes in L and to flow on toward B. With the choice A=L, we would have excluded this possibility, which would have caused a notable impact on TPT computations, given the size of our boxes. The effective currents of reactive trajectories resulting from the TPT analysis are depicted in Fig. 5, with the target set B=G indicated by the red boxes. In black, we depict the subset of pollution-capable coastal boxes L.

FIG. 5.

Inferred reactive probability currents of marine debris into garbage patches (red boxes). Black boxes indicate coastal boxes from which those currents emerge.

FIG. 5.

Inferred reactive probability currents of marine debris into garbage patches (red boxes). Black boxes indicate coastal boxes from which those currents emerge.

Close modal

We first note that the extent of the reactive currents running into the subtropical gyre patches is in general larger than that of those running into the near coastal patches. This indicates that the near coastal patches are mainly fed from nearby land-based plastic waste sources. An exception is the Bay of Bengal patch, which appears to accumulate garbage from remote sources in the coasts of the Arabian Sea and even more remote ones in the coasts of Indonesia. Particularly constrained seems to be the patch in the Gulf of Guinea, which is inferred to be filled with plastic debris releases at the southern coasts of West Africa. A refined assessment of these mostly qualitative conclusions is presented below.

Continuing with the visual inspection of Fig. 5, for the North Pacific patch, the TPT analysis infers a robust zonal eastward reactive channel into it straight out from the coasts of China. A good deal of the transported debris is inferred to travel back to the western side of the North Pacific basin, where the stationary distribution p|O also tends to maximize. A pollution source for the South Pacific patch is not restricted to the east coast of South America. Indeed, a westerly transition channel originating in the Indian Ocean and the coasts of New Zealand is also identified. Two clear reactive paths into the North Atlantic are identified, one mainly coming from the southeastern coast of the United States and another one coming from the northern coasts of West Africa. In turn, the South Atlantic patch is fed from debris transport from the Brazil–Malvinas Confluence and the southern tip of Africa. A main carrier of pollution for the Indian Ocean patch is the Agulhas Return Current. However, this pollution channel bifurcates a bit east of the patch’s longitude, where a branch originates to ultimately feed the South Pacific patch. Indeed, the pattern of the currents near the Indian Ocean patch does not suggest as clear channels into it as into the patches in the other subtropical gyres. This seems consistent with the reported55 weakness of the Indian Ocean patch.

It is important to note that while ocean currents play a dominant role in transporting debris, the reactive paths inferred by TPT do not resemble entirely the mean surface-ocean currents. However, this is not unexpected given the various mechanisms, noted above, controlling the motion of the floating material beyond advection by ocean currents. We stress again that TPT, by construction, highlights currents composed of only trajectories that go from A (source) to B (target), thereby excluding information about currents that go from B to A, A to A, and B to B.

The expected transition duration (27) is estimated to be 2.6 years from the coasts into the subtropical gyre patches and the gulf and coastal sea patches. Note that this is the mean time a reactive trajectory takes from being injected into the oceans to hit any of these patches. If we set B=Y, where YG is the set of indices corresponding to the subtropical gyre patches, then the mean duration is 5.6 years, cf. (26) and (27). The expected durations of individual transition paths into the North Pacific, South Pacific, North Atlantic, South Atlantic, and Indian Ocean patches are 7.3, 8.6, 4.3, 4.0, and 4.2 yr, respectively. The mean durations of those into the patches in the Bay of Bengal and the Gulf of Guinea are 0.6 and 0.2 years, respectively. The proximity to the coasts explains the short mean durations of the latter transition channels. As for the transition channels into the subtropical gyre patches, those into the South Pacific and North Atlantic patches stand out as the overall slowest and fastest in the class, respectively. These times to individual patches represent the mean duration of those reactive trajectories that first hit the set B=gG through the respective patch, i.e., transitions are direct and not through other patches, Gg (which can be avoided by using TPT for open dynamics).

Additional insight is provided by the normalized distribution of reactive trajectories μ^AB, as plotted in Fig. 6, showing where reactive trajectories spend most of the time while on their way from source A to target B. Note that the reactive trajectories tend to bottleneck over large regions around the subtropical gyre patches except the Indian Ocean gyre patch. These regions measure the size of the patches. The bottleneck is particularly pretty intense in the North Pacific gyre. The reactive flows in the Indian Ocean patch are not seen to spend as much time near the patch as near the other subtropical gyre patches. This is consistent with it being a weak garbage patch. Additional intense bottlenecks are observed to concentrate in the Bay of Bengal.

FIG. 6.

Probability density of reactive debris paths that indicate where debris bottlenecks on their way into the garbage patches.

FIG. 6.

Probability density of reactive debris paths that indicate where debris bottlenecks on their way into the garbage patches.

Close modal

Further insight into the domain of influence of each individual garbage patch gG, and thus into the locations on the coast where debris flowing into them originate from, is offered by associating to each state iO the most likely patch g (target) to hit according to the probability in i to forward-commit to g, viz.,

(37)

This way every box of the partition gets assigned to a patch, forming what we call a forward-committor-based dynamical geography, which is shown in Fig. 7. Note the large influence exerted by the subtropical patches on the global transport of marine debris, particularly those in the subtropical gyres whose provinces, which play a role analogous to basins of attraction, span the largest areas. Similar influence of the subtropical patches was inferred from the spectral analysis12,14 applied on simulated trajectories18 and from direct evolutions using drifter trajectory data.3 The relatively large influence of the Bay of Bengal patch inferred from the visual inspection of the reactive currents into it is well framed by the geography.

FIG. 7.

Forward-committor-based dynamical geography revealing domains of influence for the garbage patches with the provinces colored according to residence time.

FIG. 7.

Forward-committor-based dynamical geography revealing domains of influence for the garbage patches with the provinces colored according to residence time.

Close modal

The provinces of the geography in Fig. 7 are colored according to the mean residence time, defined as follows. Let QO be the box indices of a given province. The mean time it takes a trajectory initialized in iQ to move out of Q and thus hit the complement of Q, hiQ:=E(τOωQ+X0=i), is given by the solution of the linear equation,10,34,61

(38)

where hQ=(hiQ)iQ. By taking the average of hQ with respect to the stationary density p|Q, we get the residence time in Q, i.e.,

(39)

The longest residence time is 14.6 years, computed for the South Pacific province, whereas the shortest residence times are 0.7 and 0.3 years for the Bay of Bengal and Gulf of Guinea regions, respectively. The North Pacific Ocean, North Atlantic Ocean, and South Atlantic Ocean subtropical subtropical garbage patches all have comparable residence times that range between 7 and 7.5 years, while the Indian Ocean garbage patch has a much shorter residence time of 1.8 years.

The interconnectivity of the subtropical garbage patches with respect to the amount of debris particles that are exchanged between patches is presented in Fig. 8. More precisely, we compute the reactive flux from A=yYG to B=(Yy)ω, where Y is the set of subtropical gyre patches. Then, the proportions of total debris mass present in the ocean that flow per time unit (T=5d which is one time step of the Markov chain) out of A and make their way toward B are kAB=1.4×104, 3.7×105, 1.1×104, 6.2×105, and 6.3×105 for A chosen as the North Pacific, South Pacific, North Atlantic, South Atlantic, and Indian Ocean patches, respectively. We can further decompose the transition rate from A to B into the sum of arrival rates into each individual patch b in B, kAB=kB=bBkb as in (26). For a fixed A, the arrival rates into each bB are shown in the rows of Fig. 8. Consistently, the “emission” from a garbage patch y recirculates almost completely through ω before reaching any other patch Yy, hence the much higher rates in the column corresponding to the nirvana state. In addition, relatively high rates between the subtropical garbage patches of the southern hemisphere highlight an interconnection between the Indian Ocean, the South Atlantic and the South Pacific patches. Specifically, the Southern Atlantic debris transit at high rate to the South Pacific and Indian Oceans, and similarly, the Indian Ocean debris transit at high rate to the South Pacific and Atlantic Ocean garbage patches. Finally, the reactive rates from the North Atlantic gyre to any other subtropical garbage patches are negligible, confirming again that it has very little connection with other patches and debris that manage to escape it most likely end up on land or in ice.

FIG. 8.

Transition rates from each subtropical gyre patch, presented by row, into all other subtropical gyre patches and the nirvana state. The last row shows the rates from the nirvana state ω into the subtropical gyre patches from the results presented in Fig. 5.

FIG. 8.

Transition rates from each subtropical gyre patch, presented by row, into all other subtropical gyre patches and the nirvana state. The last row shows the rates from the nirvana state ω into the subtropical gyre patches from the results presented in Fig. 5.

Close modal

The last row of Fig. 8 shows transition rates from (26) corresponding to the currents into each subtropical patch presented in Fig. 5 with A chosen as the nirvana state ω. As expected, the transition rates from ω to the subtropical garbage patches are orders of magnitude higher than the transition rates between patches. Bearing in mind that those transition rates are very low, meaning that the transitions are unlikely, associated reactive currents are depicted in Figs. 9 and 10. These represent potential pathways that marine debris might take out of the gyres, for instance, in the event of unusually strong winds.

FIG. 9.

Reactive currents from the union of all subtropical gyres to the nirvana state ω. The source set is indicated in black, and the red boxes are coastal boxes with a nonzero fraction of land where debris can beach.

FIG. 9.

Reactive currents from the union of all subtropical gyres to the nirvana state ω. The source set is indicated in black, and the red boxes are coastal boxes with a nonzero fraction of land where debris can beach.

Close modal
FIG. 10.

Reactive currents from each subtropical gyre patch to all other subtropical gyre patches. The source set is indicated in black in each panel; the target sets are indicated in red. Note the difference in the scales across the panels.

FIG. 10.

Reactive currents from each subtropical gyre patch to all other subtropical gyre patches. The source set is indicated in black in each panel; the target sets are indicated in red. Note the difference in the scales across the panels.

Close modal

Figure 9 presents the reactive currents from the subtropical gyre patches to the nirvana state ω, which corresponds to the last column of Fig. 8. That is, we set A=Y (black squares) and B=ω (red squares are coastal bins iL, where POω>0). In general, debris out of the northern hemisphere patches have a larger probability of beaching than the southern hemisphere patches. In particular, the reactive currents in the Indian Ocean follow the general path of debris from the search area of the infamous Malaysia Airlines flight MH370 to the locations of recovered debris on the coasts of Mauritius, Madagascar, Mozambique, Tanzania, and South Africa.33 

In turn, Fig. 10 presents the reactive currents from a subtropical gyre patch to the union of all other subtropical gyre patches. To place the focus on debris trajectories that stay in the ocean, we do not allow reactive passages via ω. Thus, we use TPT for open domains by setting A=y (black square) and B=Yy (red squares) in (28) and (29). The reactive currents out of the Indian Ocean patch are quite strong, in agreement with reports55 on its weak character. However, these are somewhat weaker than those out of the South Atlantic patch. Note that both the Indian Ocean patch and the South Atlantic patch exchange debris with the South Pacific Ocean patch, as shown in Fig. 8, through the Antarctic Circumpolar Current. The currents that flow out of the North Pacific patch are much weaker, yet not as weak as those coming out of the North Atlantic patch. The strength of the currents out of the South Pacific patch ranges in between the above.

To quantify the above qualitative conclusions from the inspection of the transition channels, we computed the reactive rates from each yY to Yy, telling us the amount of debris probability mass that flows out of yY per time step (i.e., per five days) and is on its direct way to Yy, equal to the row sums of Fig. 8 excluding the portion that goes into nirvana. The transition rate (24) gives 5.6×109, 1.3×106, 3.7×1010, 1.4×105, and 4.5×106 for the North Pacific, South Pacific, North Atlantic, South Atlantic, and Indian Ocean patches, respectively, which confirm our qualitative assessments above. We note that each of these rates is at least one order of magnitude smaller than those reported at the very beginning of this section, except that of the South Atlantic, where it is merely a factor 5 weaker. This indicates that debris leaving the South Atlantic is most frequently finding its way to other patches.

It must be noted that the above transition rates do not say anything about retention. They tell us which patch “emits” the most frequently such debris that finds its way to another patch. A low rate does not need to mean that debris leaving A come back to A since the debris can hit ω too before hitting B, as shown by the much higher transition rates to the nirvana state in Fig. 8. In other words, a low rate should not be taken to mean the same as high attraction. Thus, the transition rate computation results just described do not contradict those from direct density evolution in Fig. 4, which had identified the North Pacific patch as the most attracting of all.

We have presented a novel application of transition path theory (TPT), here extended to open autonomous dynamical systems. The problem chosen was that of pollution routes from possible coastline sources (reactive states) into garbage patches in the global surface ocean (product states).

Undrogued drifter trajectories from the NOAA Global Drifter Program were used to derive a Markov chain on which TPT was applied, as a model for the time-asymptotic dynamics of marine debris pollution. Modeling the probability of trajectories to beach as a function of the fraction of land filling each coastal box of the covering of the world ocean domain resulted in an open system, which was closed by sending the probability imbalance back into the chain according to the capacity of coastal boxes to “pollute the oceans” as measured by its share of global mismanaged plastic waste. Assuming a constant pollution rate, our time-homogeneous model was shown to be the statistical equivalent of a “saturated” (stationary) pollution redistribution dynamics.

A high probability transition channel was identified connecting the Great Pacific garbage patch with the coasts of Eastern Asia, suggesting an important source of plastic pollution there. The weakness of the Indian Ocean gyre as a trap of plastic debris was found consistent with transition paths not converging in the gyre. While the North Pacific subtropical gyre was found to be most attracting consistent with earlier assessments, the South Pacific gyre stood out as the most enduring in the sense that the total reactive rate out of that gyre into other gyres and the nirvana state resulted the smallest of all. The weakest of all the gyres in terms of its capacity to trap and hold within plastic waste resulted to be the South Atlantic gyre. The gyres were found in general weakly communicated. Indeed, in the event of anomalously intense winds, a subtropical gyre is more likely to export garbage out toward the coastlines than into another gyre.

Our results, including prospects for garbage patches yet to be directly and/or robustly observed, namely, the Gulf of Guinea and the Bay of Bengal, have implications for activities such as ocean cleanup as the revealed reactive pollution routes provide targets, alternative to the great garbage patches themselves, to aim those efforts. Additional ocean applications of TPT are under way (e.g., using submerged float data and targeting meridional overturning routes) and will be reported elsewhere.

See the supplementary material for versions of Figs. 4 and 5 assuming α=12 (Figs. S1 and S4, respectively), 34 (Figs. S2 and S5), and 1 (Figs. S3 and S6) in (32).

We thank María J. Olascoaga for the benefit of many discussions on Lagrangian ocean dynamics. Support for this work was provided by National Science Foundation (NSF) Grant No. OCE1851097 (P.M. and F.J.B.-V.), Deutsche Forschungsgemeinschaft (DFG) through Grant No. CRC 1114 “Scaling Cascades in Complex Systems,” Project No. 235221301, Project A01 (P.K.), and Germany’s Excellence Strategy—The Berlin Mathematics Research Center MATH+ (EXC-2046/1, Project ID: 390685689) (L.H.).

The data that support the findings of this study are openly available in the NOAA Global Drifter Program dataset at http://www.aoml.noaa.gov/phod/dac/, Ref. 62. The mismanaged plastic waste data are distributed from https://ourworldindata.org/, Ref. 53. Finally, the numerical code to reproduce the findings of this study is openly available at https://github.com/philippemiron/pygtm, Ref. 63.

1.
R.
Lumpkin
and
M.
Pazos
, “Measuring surface currents with surface velocity program drifters: The instrument, its data and some recent results,” in Lagrangian Analysis and Prediction of Coastal and Ocean Dynamics, edited by A. Griffa, A. D. Kirwan, A. Mariano, T. Özgökmen, and T. Rossby (Cambridge University Press, 2007), Chap. 2, pp. 39–67.
2.
E.
van Sebille
,
E. H.
England
, and
G.
Froyland
, “
Origin, dynamics and evolution of ocean garbage patches from observed surface drifters
,”
Environ. Res. Lett.
7
,
044040
(
2012
).
3.
A. N.
Maximenko
,
J.
Hafner
, and
P.
Niiler
, “
Pathways of marine debris derived from trajectories of Lagrangian drifters
,”
Mar. Pollut. Bull.
65
,
51
62
(
2012
).
4.
A.
Cozar
,
F.
Echevarria
,
J. I.
Gonzalez-Gordillo
,
X.
Irigoien
,
B.
Ubeda
,
S.
Hernandez-Leon
,
A. T.
Palma
,
S.
Navarro
,
J.
Garcia-de Lomas
,
R.
andrea
,
M. L.
Fernandez-de Puelles
, and
C. M.
Duarte
, “
Plastic debris in the open ocean
,”
Proc. Natl. Acad. Sci. U.S.A.
111
,
10239
10244
(
2014
).
5.
F. J.
Beron-Vera
,
M. J.
Olascoaga
, and
R.
Lumpkin
, “
Inertia-induced accumulation of flotsam in the subtropical gyres
,”
Geophys. Res. Lett.
43
,
12228
12233
, https://doi.org/10.1002/2016GL071443 (
2016
).
6.
F. J.
Beron-Vera
,
M. J.
Olascoaga
, and
P.
Miron
, “
Building a Maxey–Riley framework for surface ocean inertial particle dynamics
,”
Phys. Fluids
31
,
096602
(
2019
).
7.
F. J.
Beron-Vera
, “Nonlinear dynamics of inertial particles in the ocean: From drifters and floats to marine debris and Sargassum,” Nonlinear Dyn. 103, 1–26 (2021). doi:10.1007/s11071-020-06053-z.
8.
A.
Lasota
and
M. C.
Mackey
, Chaos, Fractals and Noise: Stochastic Aspects of Dynamics, 2nd ed., Applied Mathematical Sciences Vol. 97 (Springer, New York, 1994).
9.
P.
Brémaud
, Markov Chains, Gibbs Fields Monte Carlo Simulation Queues, Texts in Applied Mathematics, Vol. 31 (Springer, New York, 1999).
10.
J.
Norris
,
Markov Chains
(
Cambridge University Press
,
1998
).
11.
M.
Dellnitz
and
A.
Hohmann
, “
A subdivision algorithm for the computation of unstable manifolds and global attractors
,”
Numer. Math.
75
,
293
317
(
1997
).
12.
M.
Dellnitz
and
O.
Junge
, “
On the approximation of complicated dynamical behavior
,”
SIAM J. Numer. Anal.
36
,
491
515
(
1999
).
13.
G.
Froyland
and
M.
Dellnitz
, “
Detecting and locating near-optimal almost-invariant sets and cycles
,”
SIAM J. Sci. Comput.
24
,
1839
1863
(
2003
).
14.
P.
Koltai
, “Efficient approximation methods for the global long-term behavior of dynamical systems—Theory, algorithms and examples,” Ph.D. thesis (Technical University of Munich, 2010).
15.
E.
Ser-Giacomi
,
V.
Rossi
,
C.
Lopez
, and
E.
Hernandez-Garcia
, “
Flow networks: A characterization of geophysical fluid transport
,”
Chaos
25
,
036404
(
2015
).
16.
M.
Dubois
,
V.
Rossi
,
E.
Ser-Giacomi
,
S.
Arnaud-Haond
,
C.
López
, and
E.
Hernández-García
, “
Linking basin-scale connectivity, oceanography and population dynamics for the conservation and management of marine ecosystems
,”
Glob. Ecol. Biogeogr.
25
,
503
515
(
2016
).
17.
P.
Monroy
,
E.
Hernández-García
,
V.
Rossi
, and
C.
López
, “
Modeling the dynamical sinking of biogenic particles in oceanic flow
,”
Nonlinear Process. Geophys.
24
,
293
305
(
2017
).
18.
G.
Froyland
,
R. M.
Stuart
, and
E.
van Sebille
, “
How well-connected is the surface of the global ocean?
,”
Chaos
24
,
033126
(
2014
).
19.
E.
Vanden-Eijnden
, “
Transition path theory
,” in
Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1. Lecture Notes in Physics
,
edited by M. Ferrario, G. Ciccotti, and K. Binder (Springer, Berlin, Heidelberg, 2006), Vol. 703
.
20.
P.
Metzner
,
C.
Schütte
, and
E.
Vanden-Eijnden
, “
Illustration of transition path theory on a collection of simple examples
,”
J. Chem. Phys.
125
,
084110
(
2006
).
21.
E.
Weinan
and
E.
Vanden-Eijnden
, “
Towards a theory of transition paths
,”
J. Stat. Phys.
123
,
503
623
(
2006
).
22.
P.
Metzner
,
C.
Schütte
, and
E.
Vanden-Eijnden
, “
Transition path theory for Markov jump processes
,”
Multiscale Model. Simul.
7
,
1192
1219
(
2009
).
23.
E.
Weinan
and
E.
Vanden-Eijnden
, “
Transition-path theory and path-finding algorithms for the study of rare events
,”
Annu. Rev. Phys. Chem.
61
,
391
420
(
2010
).
24.
E.
Wigner
, “
The transition state method
,”
Trans. Faraday Soc.
34
,
29
41
(
1938
).
25.
L.
Pratt
, “
A statistical method for identifying transition states in high dimensional problems
,”
J. Chem. Phys.
85
,
5045
5048
(
1986
).
26.
F.
Noé
,
C.
Schütte
,
E.
Vanden-Eijnden
,
L.
Reich
, and
T. R.
Weikl
, “
Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations
,”
Proc. Natl. Acad. Sci. U.S.A.
106
,
19011
19016
(
2009
).
27.
V. A.
Voelz
,
G. R.
Bowman
,
K.
Beauchamp
, and
V. S.
Pande
, “
Molecular simulation of ab initio protein folding for a millisecond folder ntl9 (1- 39)
,”
J. Am. Chem. Soc.
132
,
1526
1528
(
2010
).
28.
D.
Lucente
,
S.
Duffner
,
C.
Herbert
,
J.
Rolland
, and
F.
Bouchet
, “Machine learning of committor functions for predicting high impact climate events,” arXiv:1910.11736 (2019).
29.
J.
Finkel
,
D. S.
Abbot
, and
J.
Weare
, “
Path properties of atmospheric transitions: Illustration with a low-order sudden stratospheric warming model
,”
J. Atmos. Sci.
77
,
2327
2347
(
2020
).
30.
G.
Froyland
,
P. K.
Pollett
, and
R. M.
Stuart
, “
A closing scheme for finding almost-invariant sets in open dynamical systems
,”
J. Comput. Dyn.
1
,
135
(
2014
).
31.
B.
Lünsmann
and
H.
Kantz
, “
An extended transfer operator approach to identify separatrices in open flows
,”
Chaos
28
,
053101
(
2018
).
32.
P.
Miron
,
F. J.
Beron-Vera
,
M. J.
Olascoaga
,
J.
Sheinbaum
,
P.
Pérez-Brunius
, and
G.
Froyland
, “
Lagrangian dynamical geography of the Gulf of Mexico
,”
Sci. Rep.
7
,
7021
(
2017
).
33.
P.
Miron
,
F. J.
Beron-Vera
,
M. J.
Olascoaga
, and
P.
Koltai
, “
Markov-chain-inspired search for MH370
,”
Chaos
29
,
041105
(
2019
).
34.
P.
Miron
,
F. J.
Beron-Vera
,
M. J.
Olascoaga
,
G.
Froyland
,
P.
Pérez-Brunius
, and
J.
Sheinbaum
, “
Lagrangian geography of the deep gulf of Mexico
,”
J. Phys. Oceanogr.
49
,
269
290
(
2019
).
35.
M. J.
Olascoaga
,
P.
Miron
,
C.
Paris
,
P.
Pérez-Brunius
,
R.
Pérez-Portela
,
R. H.
Smith
, and
A.
Vaz
, “
Connectivity of pulley ridge with remote locations as inferred from satellite-tracked drifter trajectories
,”
J. Geophys. Res.
123
,
5742
5750
, https://doi.org/10.1029/2018JC014057 (
2018
).
36.
F. J.
Beron-Vera
,
N.
Bodnariuk
,
M.
Saraceno
,
M. J.
Olascoaga
, and
C.
Simionato
, “
Stability of the Malvinas current
,”
Chaos
30
,
013152
(
2020
).
37.
E.
Morrison
,
A.
Shipman
,
S.
Shrestha
,
E.
Squier
, and
K.
Stack Whitney
, “
Evaluating the ocean cleanup, a marine debris removal project in the North Pacific Gyre, using SWOT analysis
,”
Case Stud. Environ.
3
,
1
6
(
2019
).
38.
S. M.
Ulam
, A Collection of Mathematical Problems, Interscience Tracts in Pure and Applied Mathematics (Interscience, 1960).
39.
Z.
Kovács
and
T.
Tél
, “
Scaling in multifractals: Discretization of an eigenvalue problem
,”
Phys. Rev. A
40
,
4641
4646
(
1989
).
40.
L.
Helfmann
,
E. R.
Borrell
,
C.
Schütte
, and
P.
Kotai
, “
Extending transition path theory: Periodically driven and finite-time dynamics
,”
J. Nonlinear Sci.
30
,
3321
3366
(
2020
).
41.
A. L.
Sybrandy
and
P. P.
Niiler
, “WOCE/TOGA Lagrangian drifter contruction manual,” Technical Report SIO Reference 91/6 (Scripps Institution of Oceanography, La Jolla, CA, 1991).
42.
P. P.
Niiler
and
J. D.
Paduan
, “
Wind-driven motions in the northeastern Pacific as measured by Lagrangian drifters
,”
J. Phys. Oceanogr.
25
,
2819
2830
(
1995
).
43.
R.
Lumpkin
,
S. A.
Grodsky
,
L.
Centurioni
,
M.-H.
Rio
,
J. A.
Carton
, and
D.
Lee
, “
Removing spurious low-frequency variability in drifter velocities
,”
J. Atmos. Oceanic. Technol.
30
,
353
360
(
2012
).
44.
M. J.
Olascoaga
,
F. J.
Beron-Vera
,
P.
Miron
,
J.
Triñanes
,
N. F.
Putman
,
R.
Lumpkin
, and
G. J.
Goni
, “
Observation and quantification of inertial effects on the drift of floating objects at the ocean surface
,”
Phys. Fluids
32
,
026601
(
2020
).
45.
P.
Miron
,
S.
Medina
,
M. J.
Olascaoaga
, and
F. J.
Beron-Vera
, “
Laboratory verification of a Maxey–Riley theory for inertial ocean dynamics
,”
Phys. Fluids
32
,
071703
(
2020
).
46.
M.
Sarich
,
R.
Banisch
,
C.
Hartmann
, and
C.
Schütte
, “
Markov state models for rare events in molecular dynamics
,”
Entropy
16
,
258
286
(
2014
).
47.
R.
McAdam
and
E.
van Sebille
, “
Surface connectivity and interocean exchanges from drifter-based transition matrices
,”
J. Geophys. Res.: Oceans
123
,
514
532
, https://doi.org/10.1002/2017JC013363 (
2018
).
48.
A.
Gordon
and
R.
Fine
, “
Pathways of water between the Pacific and Indian oceans in the Indonesian seas
,”
Nature
379
,
46
149
(
1996
).
49.
D.
Tillinger
and
A. L.
Gordon
, “
Fifty years of the Indonesian throughflow
,”
J. Clim.
22
,
6342
6355
(
2009
).
50.
R.
Tarjan
, “
Depth-first search and linear graph algorithms
,”
SIAM J. Comput.
1
,
146
160
(
1972
).
51.
M.
Kubota
, “
A mechanism for the accumulation of floating marine debris north of Hawaii
,”
J. Phys. Oceanogr.
24
,
1059
1064
(
1994
).
52.
L.
Lebreton
,
B.
Slat
,
F.
Ferrari
,
B.
Sainte-Rose
,
J.
Aitken
,
R.
Marthouse
,
S.
Hajbane
,
S.
Cunsolo
,
A.
Schwarz
,
A.
Levivier
,
K.
Noble
,
P.
Debeljak
,
H.
Maral
,
R.
Schoeneich-Argent
,
R.
Brambini
, and
J.
Reisser
, “
Evidence that the great pacific garbage patch is rapidly accumulating plastic
,”
Sci. Rep.
8
,
4666
(
2018
).
53.
J. R.
Jambeck
,
R.
Geyer
,
C.
Wilcox
,
T. R.
Siegler
,
M.
Perryman
,
A.
Andrady
,
R.
Narayan
, and
K. L.
Law
, “
Plastic waste inputs from land into the ocean
,”
Science
347
,
768
771
(
2015
).
54.
R.
Yamashita
and
A.
Tanimura
, “
Floating plastic in the Kuroshio Current area, western North Pacific Ocean
,”
Mar. Pollut. Bull.
54
,
485
488
(
2007
).
55.
M.
van der Mheen
,
C.
Pattiaratchi
, and
E.
van Sebille
, “
Role of indian ocean dynamics on accumulation of buoyant debris
,”
J. Geophys. Res.: Oceans
124
,
2571
2590
, https://doi.org/10.1029/2018JC014806 (
2019
).
56.
P.
Miron
,
M. J.
Olascoaga
,
F. J.
Beron-Vera
,
J.
Triñanes
,
N. F.
Putman
,
R.
Lumpkin
, and
G. J.
Goni
, “
Clustering of marine-debris-and Sargassum-like drifters explained by inertial particle dynamics
,”
Geophys. Res. Lett.
47
,
e2020GL089874
, https://doi.org/10.1029/2020GL089874 (
2020
).
57.
A.
Cozar
,
E.
Marti
,
C. M.
Duarte
,
J.
Garcia-de Lomas
,
E.
van Sebille
,
T. J.
Ballatore
,
V. M.
Eguiluz
,
J. I.
Gonzalez-Gordillo
,
M. L.
Pedrotti
,
F.
Echevarria
,
R.
Trouble
, and
X.
Irigoien
, “
The arctic ocean as a dead end for floating plastics in the north atlantic branch of the thermohaline circulation
,”
Sci. Adv.
3
,
e1600582
(
2017
).
58.
A.
Mountford
and
M. A.
Morales Maqueda
, “
Eulerian modelling of the three-dimensional distribution of seven popular microplastic types in the global ocean
,”
J. Geophys. Res.: Oceans
124
,
8558
8573
, https://doi.org/10.1029/2019JC015050 (
2019
).
59.
T. T.
Sutton
,
M. R.
Clark
,
D. C.
Dunn
,
P. N.
Halpin
,
A. D.
Rogers
,
J.
Guinotte
,
S. J.
Bograd
,
M. V.
Angel
,
J. A. A.
Perez
,
K.
Wishner
,
R. L.
Haedrich
,
D. J.
Lindsay
,
J. C.
Drazen
,
A.
Vereshchaka
,
U.
Piatkowski
,
T.
Morato
,
K.
Blachowiak-Samolyk
,
B. H.
Robison
,
K. M.
Gjerde
,
A.
Pierrot-Bults
,
P.
Bernal
,
G.
Reygondeau
, and
M.
Heino
, “
A global biogeographic classification of the mesopelagic zone
,”
Deep Sea Res.
126
,
85
102
(
2017
).
60.
P.
Ryan
, “
A simple technique for counting marine debris at sea reveals steep litter gradients between the Straits of Malacca and the Bay of Bengal
,”
Mar. Pollut. Bull.
69
,
128
136
(
2013
).
61.
M.
Dellnitz
,
G.
Froyland
,
C.
Horenkam
,
K.
Padberg-Gehle
, and
A.
Sen Gupta
, “
Seasonal variability of the subpolar gyres in the southern ocean: A numerical investigation based on transfer operators
,”
Nonlinear Process. Geophys.
16
,
655
663
(
2009
).
62.
R.
Lumpkin
, and
L. Centurioni
Luca
, “
Global Drifter Program quality-controlled 6-hour interpolated data from ocean surface drifting buoys
,” NOAA National Centers for Environmental Information (2019).
Dataset.
Accessed Oct. 1, 2021. (
2019
).
63.
P.
Miron
, “
python Geospatial Transition Matrix toolbox (pyGTM)
,” Accessed on October (
2020
). https://github.com/philippemiron/pygtm.

Supplementary Material