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.

## I. INTRODUCTION

The long-term fate of satellite-tracked drifting buoys from the National Oceanic and Atmospheric Administration's (NOAA) Global Drifter Program^{1} is characterized by a tendency to form clusters in the oceans’ subtropical gyres^{2,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 explained^{5–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 theory^{8} and Markov chains,^{9,10} which form the basis for approximating asymptotically invariant sets using so-called *set-oriented methods*^{11–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 theory^{24} and transition path sampling.^{25} Since then, the TPT framework has also been applied to study molecular conformation changes^{26,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 approaches^{30,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.

## II. SETUP FOR CLOSED DYNAMICAL SYSTEMS

Let us assume that floating debris trajectories are described by a time-homogeneous stochastic process in continuous space $X\u2282R2$ and observed at discrete times $nT$, $n\u2208Z$. Its transition probabilities are controlled by a stochastic kernel $K(x,y)\u22650$ such that $\u222bXK(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 $\sigma $-algebra of subsets measured by (normalized) area. Then, a probability density $f(x)\u22650$, $\u222bXf(x)dx=1$, describing the distribution of the random position $XnT$ at any time $nT$ evolves to the distribution

at time $(n+1)T$, which defines a Markov operator $P:L1(X)\u21ba$ 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}i\u2208S$, $S:={1,\u2026,N}\u2282Z+$, 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)}i\u2208S$, where $1A(x)=1$ if $x\u2208A$ and 0 otherwise. The discrete action of $P$ on $VN$ is described by a matrix $P=(Pij)i,j\u2208S\u2208RN\xd7N$ called a *transition matrix*. The transition matrix results from the projection,^{33,34}

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.,

Note that $\u2211j\u2208SPij=1$ for all $i\u2208S$, 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)i\u2208S$, $\u2211i\u2208Sfi=1$, is calculated under *left* multiplication, i.e.,

as it follows by noting that $Pr(X(n+1)T\u2208Bj)=\u2211i\u2208SPr(X(n+1)T\u2208Bj,XnT\u2208Bi)=\u2211i\u2208SPr(XnT\u2208Bi)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,\u2026,1)$ is a *right* eigenvector with eigenvalue $\lambda =1$, i.e., $P1=1$. The eigenvalue $\lambda =1$ is the largest eigenvalue of $P$. The associated potentially nonunique *left* eigenvector $p=(pi)i\u2208S$ is invariant, because $pP=p$ and can be chosen componentwise nonnegative (by the Perron–Frobenius theorem).

We call $P$ *irreducible* (or *ergodic*) if for all $i,j\u2208S$ there exists $nij\u2208Z0+\u2216{\u221e}$ such that $(Pnij)ij>0$. To wit, all states of an irreducible Markov chain communicate, the eigenvalue $\lambda =1$ is simple, and the corresponding left eigenvector $p$ is strictly positive.^{10} We call $P$ *aperiodic* (or *mixing*) if there exists $i\u2208S$ such that $gcd{n\u2208Z0+:(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 ($\u2211i\u2208Spi=1$), satisfies $0<p=pP=limn\u2191\u221efPn$ for any probability vector $f$. We call $p$ an invariant limiting probability vector or *stationary distribution*.

We adopt the traditional notation with ${Xt}t\u2208Z$ instead of ${XnT}n\u2208Z$ and write, for instance, $Pij=Pr(Xt+1=j\u2223Xt=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(Xt\u2208Bi)=pi$ for all $t\u2208Z$.

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.

## III. 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}i\u2208O$ of the ocean domain $X$ but the probability to transition from one box with index $i\u2208O$ 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=j\u2223Xt=i)$ for $i,j\u2208O$. Since the rows of $PO$ no longer have to add up to one, $PO$ represents a *substochastic* matrix.

We assume that a larger domain $S\u2283O$ exists on which the dynamics are closed, i.e., the transition matrix $P$ on box entries $i,j\u2208S$ 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)i\u2208S$, while we denote the restriction to the open domain by $p|O=(pi)i\u2208O$.

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 $\omega $, which we will call a *two-way nirvana state*, and letting all the outflow from $O$ flow into $\omega $, while also redistributing the probability mass from $\omega $ *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 ${\omega}$ also by $\omega $ and refer to it too as the two-way nirvana state.

The resulting transition matrix on $O\u222a\omega $ reads (possibly overloading the notation by denoting it by $P$ again)

where $PO\u2192\omega :=(1\u2212\u2211j\u2208OPijO)i\u2208O$ (understood as a column vector) gives the outflow from $O$ to $\omega $ and $P\omega \u2192O$ is a (row) vector that gives the inflow and has to be a probability vector. Note that the matrix $P$ is stochastic $\u2211j\u2208O\u222a\omega Pij=1$ for all $i\u2208O\u222a\omega $ 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 choice^{30} for $P\omega \u2192O$ is to redistribute according to the quasistationary distribution of $PO$. Lünsmann and Kantz^{31} 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 $O\u228aS$. 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.

## IV. TRANSITION PATH THEORY

### A. TPT for closed systems

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 $A\u2282S$ to another, disjoint set $B\u2282S$ 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).

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 $S\u2282S$ is the stopping time random variable defined as

where $inf\u2205:=\u221e$. The *forward committor*$q+:=(qi+)i\u2208S$ gives the probability that a trajectory starting in $i\u2208S$*first enters*$B$, *not*$A$, i.e.,

Note that $qi\u2208A+=0$ while $qi\u2208B+=1$. For $i\u2208C:=S\u2216(A\u222aB)$, one has that

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

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

The *last exit time*, in turn, is defined by

where $sup\u2205:=\u2212\u221e$, which is a stopping time, but for the *time-reversed* chain ${Xt\u2212}t\u2208Z$ that traverses the original Markov chain backward in time, i.e., $Xt\u2212:=X\u2212t$. The reversed chain’s transition matrix, $P\u2212=(Pij\u2212)i,j\u2208S$, is given by

since the chain is assumed to be in stationarity. The time-reversed transition matrix $P\u2212$ is ergodic and mixing and has the same stationary distribution $p$ as $P$. The *backward committor*$q\u2212:=(qi\u2212)i\u2208S$ gives the probability that a trajectory starting in $i\u2208S$*last exits*$A$, *not*$B$,

In this case,

for $i\u2208C$, subject to $qi\u2208B\u2212=0$ and $qi\u2208A\u2212=1$. The (unique) solution in matrix notation,

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

The committors contain information that enable the computation of various transition statistics. The *distribution of reactive trajectories* $\mu AB=(\mu iAB)i\u2208S$, defined as the joint probability that the chain is in state $i$ while *transitioning from* $A$ *to* $B$, viz.,

tells us where reactive trajectories spend most of their time. Note that $\mu i\u2208A\u222aBAB=0$. The distribution of reactive trajectories is computable from the committor probabilities and the stationary distribution,

A *density of reactive trajectories* $\mu ^AB=(\mu ^iAB)i\u2208S$ is obtained by normalizing $\mu iAB$ by the probability to be reactive,

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

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 trajectories* $fAB=(fijAB)i,j\u2208S$ gives the average flux of trajectories going through $i$ and $j$ at two consecutive times while on their way from $A$ to $B$,

which is computable as

Note that the reactive current can include direct transitions from $i\u2208A$ to $j\u2208B$, 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 trajectories* $f+=(fij+)i,j\u2208S$, which gives the net amount of reactive current going through $i$ and $j$ consecutively, viz.,

To visualize $f+$ on a flow domain covered by boxes ${Bi}i\u2208S$, 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 $\u2211j\u2260ifij+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.,

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$,

By a simple calculation, it can be shown that summing the reactive current out of $A$, $\u2211i\u2208A,j\u2208SfijAB$, is equal to aggregating the reactive current into $B$, $\u2211i\u2208S,j\u2208BfijAB$, thus

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)\u22121$-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

into the individual arrival rates into disjoint subsets $Bn$ that together give $B=\u222anBn$,

The same can also be done for decomposing $kA\u2192$.

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

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.

### B. TPT for open domains

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 $\omega $ to $O$ closes the system artificially (as in Sec. III), but we are still only interested in the transitions from $A\u2282O$ to $B\u2282O$ that *stay* in $O$ during the transition. Thus, the reactive trajectories we consider go from $A$ to $B$ without passing $A$, $B$, or $\omega $ 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 $\omega $.

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.,

while the backward committor gives the probability to have last come from $A$, not $B\u222a\omega $,

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 $\omega $ in between.

By definition, the forward committor is $qi+=0$ for $i\u2208A\u222a\omega $ and $1$ for $i\u2208B$, while in the transition region $C:=O\u2216(A\u222aB)$, it satisfies

since $q\omega +=0$ and $P$ on entries of $O$ reduces to $PO$.

The backward committor $qi\u2212=0$ for $i\u2208B$ and 1 for $i\u2208A\u222a\omega $, while, by a similar reasoning as above, it satisfies

for $i\u2208C$, where $PO,\u2212$ is the restriction of the backward-in-time transition matrix $P\u2212$ to $O$ and has entries $PijO,\u2212=pjpiPjiO$ for $i,j\u2208O$.

Therefore, system (9) remains the same with the replacement of $P$ with $PO$ and $A$ with $A\u222a\omega $. In turn, system (14) remains the same with the replacement of $P\u2212$ with $PO,\u2212$ and $B$ with $B\u222a\omega $.

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 $\mu iAB=0$ for $i=\omega $ and $fijAB=0$ for $i,j=\omega $. 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$.

## V. MARKOV-CHAIN MODEL FOR OCEAN POLLUTION

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.

### A. Preparation of *P* from drifter trajectory data

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}i\u2208S$ by lying down on the world ocean domain a grid of roughly $3\xb0$ 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 $\lambda (P(nT))\u2248\lambda (P(T))n$ is passed well up to $n=4$, where $\lambda (\u22c5)$ 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.

*Zonal connectivity*. This is addressed by identifying and continuating trajectories crossing the antimeridian connecting the eastern and western hemispheres.*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 earlier^{47}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 $j\u2208S$ such that $Bj$ is either in the Pacific Ocean or in the Atlantic Ocean.*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.*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 algorithm^{50}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.

### B. Pollution-aware model derivation

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 $\u2113:O\u2192[0,1)$ be a *land fraction* function giving the ratio between the land area and the total box area. Namely, $0<\u2113(i)<1$ for $i\u2208L\u2282O$ corresponding to boxes filled with some portion of land (or ice) (Fig. 2, top panel) and $\u2113(i)=0$ otherwise. We then follow Miron *et al.*^{33} and replace $PijO$ with

for all $i,j\u2208O$. 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 $\alpha $, 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 $\alpha =14$ produces results most consistent with them. If $\alpha =1$ (as in Miron *et al.*,^{33}) then the so-called Great Pacific garbage patch^{51} 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 $\alpha $-value assumed, as we show in the supplementary material.

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

Here, $Wi$ is the mass of mismanaged plastic waste in $Bi$, $i\u2208L$, as inferred from estimates^{53} 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 $\omega $ 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.

### C. Physical interpretation of the model

We model the distribution of garbage input per time unit to the oceans by a time-independent vector $P\omega \u2192O=:r\u2208R|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 $i\u2209L$, and $r$ is a probability vector, i.e., $\u2211iri=1$.

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

Recall that $PO$ is assumed to be irreducible; thus, $Id\u2212PO$ 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\u2192\omega $ being the vector of absorption probabilities from the boxes into nirvana (the outflow from the open domain), the stationary distribution of our chain satisfies

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

Since $\rho $ is scalar, this readily means that $p|O\u221dr(Id\u2212PO)\u22121$, 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.

## VI. LONG-TIME ASYMPTOTICS

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\omega $, 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.

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 analyses^{2,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 evidence^{4,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\omega $ 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 patch^{4} attributed to unique oceanic and atmospheric dynamics in the region^{55} 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\omega $ 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 effects^{3} 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.

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 study^{55} 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.

## VII. REACTIVE DEBRIS PATHS

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.

### A. Pollution paths into garbage patches

To infer the pollution paths into garbage patches, we choose the nirvana state $\omega $ 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=\omega $ 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$.

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 reported^{55} 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 $Y\u2282G$ 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=g\u2208G$ through the respective patch, i.e., transitions are direct and not through other patches, $G\u2216g$ (which can be avoided by using TPT for open dynamics).

Additional insight is provided by the normalized distribution of reactive trajectories $\mu ^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.

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

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 analysis^{12,14} applied on simulated trajectories^{18} 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.

The provinces of the geography in Fig. 7 are colored according to the mean *residence time*, defined as follows. Let $Q\u2282O$ be the box indices of a given province. The mean time it takes a trajectory initialized in $i\u2208Q$ to move out of $Q$ and thus hit the complement of $Q$, $hiQ:=E(\tau O\u222a\omega \u2216Q+\u2223X0=i)$, is given by the solution of the linear equation,^{10,34,61}

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

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.

### B. Pollution paths out of subtropical garbage patches

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=y\u2208Y\u2282G$ to $B=(Y\u2216y)\u222a\omega $, 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\xd710\u22124$, $3.7\xd710\u22125$, $1.1\xd710\u22124$, $6.2\xd710\u22125$, and $6.3\xd710\u22125$ 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\u2190=\u2211b\u2208Bkb\u2190$ as in (26). For a fixed $A$, the arrival rates into each $b\u2208B$ are shown in the rows of Fig. 8. Consistently, the “emission” from a garbage patch $y$ recirculates almost completely through $\omega $ before reaching any other patch $Y\u2216y$, 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.

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 $\omega $. As expected, the transition rates from $\omega $ 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.

Figure 9 presents the reactive currents from the subtropical gyre patches to the nirvana state $\omega $, which corresponds to the last column of Fig. 8. That is, we set $A=Y$ (black squares) and $B=\omega $ (red squares are coastal bins $i\u2208L$, where $PO\u2192\omega >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 $\omega $. Thus, we use TPT for open domains by setting $A=y$ (black square) and $B=Y\u2216y$ (red squares) in (28) and (29). The reactive currents out of the Indian Ocean patch are quite strong, in agreement with reports^{55} 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 $y\u2208Y$ to $Y\u2216y$, telling us the amount of debris probability mass that flows out of $y\u2208Y$ per time step (i.e., per five days) and is on its direct way to $Y\u2216y$, equal to the row sums of Fig. 8 excluding the portion that goes into nirvana. The transition rate (24) gives $5.6\xd710\u22129$, $1.3\xd710\u22126$, $3.7\xd710\u221210$, $1.4\xd710\u22125$, and $4.5\xd710\u22126$ 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 $\omega $ 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.

## VIII. SUMMARY AND CONCLUSIONS

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.).

## DATA AVAILABILITY

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.

## REFERENCES

*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.

**103**, 1–26 (2021). doi:10.1007/s11071-020-06053-z.

*Chaos, Fractals and Noise: Stochastic Aspects of Dynamics*, 2nd ed., Applied Mathematical Sciences Vol. 97 (Springer, New York, 1994).

*Markov Chains, Gibbs Fields Monte Carlo Simulation Queues*, Texts in Applied Mathematics, Vol. 31 (Springer, New York, 1999).

*A Collection of Mathematical Problems*, Interscience Tracts in Pure and Applied Mathematics (Interscience, 1960).