Transition paths of marine debris and the stability of the garbage patches
TTransition paths of marine debris and the stability of the garbage patches
P. Miron, a) F. J. Beron-Vera, L. Helfmann,
2, 3, 4 and P. Koltai Department of Atmospheric Sciences, Rosenstiel School of Marine and Atmospheric Science, University of Miami,Miami, Florida, USA Institute of Mathematics, Freie Universität Berlin, Berlin, Germany Department of Modeling and Simulation of Complex Processes, Zuse-Institute Berlin, Berlin,Germany Complexity Science Department, Potsdam Institute for Climate Impact Research, Potsdam,Germany (Dated: 24 September 2020)
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 fromtrajectories produced by satellite-tracked undrogued buoys from the NOAA Global Drifter Program. Thelatter involved coping with the openness of the system in physical space, which further required an adaptationof the standard TPT setting. Directly connecting pollution sources along coastlines with garbage patches ofvaried 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 PacificGarbage Patch; characterizing the weakness of the Indian Ocean gyre as a trap for plastic waste; and unveilinga tendency of the subtropical gyres to export garbage toward the coastlines rather than to other gyres in theevent of anomalously intense winds.PACS numbers: 02.50.Ga; 47.27.De; 92.10.Fj
Given a Markov chain, namely, a model de-scribing the stochastic state transitions in whichthe transition probability of each state dependsonly on the state attained in the previous event,transition path theory (TPT) provides a rigor-ous approach to study the statistics of transi-tions from a set of states to another, possiblydisconnected set of states. Envisioning the mo-tion of floating debris as described by a Markovchain that accounts for the ability of coastalstates to “pollute the oceans,” TPT is employedto unveil “reactive” pathways representing directtransitions from potential release locations alongthe shorelines to accumulation sites across theworld ocean. These include the subtropical gyres,whose strength in this context is investigated.
I. INTRODUCTION
The long-term fate of satellite-tracked drifting buoysfrom the NOAA Global Drifter Program is characterizedby a tendency to form clusters in the oceans’ subtropicalgyres that resemble great garbage patches. The de-velopment of such clusters, most evidently in the case ofundrogued (i.e., without a sea anchor) drifters, has beenexplained as the result of the combined action on thedrifters of converging ocean currents and winds mediatedby their inertia, which prevent them from adapting theirvelocities to that of the carrying water–air flow system. a) Electronic mail: [email protected]
The tendency of the drifters to cluster in the long runenables a probabilistic description of their dynamics us-ing results from ergodic theory and Markov chains, which form the basis for approximating asymptoticallyinvariant sets using so-called set-oriented methods . This approach places the focus on the evolution of proba-bility densities, which, unlike individual trajectories, rep-resent robust features of the dynamics. Central to thismeasure-theoretic characterization is the transfer oper-ator and the transition matrix, its discrete version re-sulting by covering the phase space with boxes, whichrepresent the states of the associated Markov chain.Such a probabilistic description has been applied on simulated drifter trajectories, suggesting a characteri-zation of great garbage patches as almost-invariant at-tracting sets with corresponding basins of attractionspanning areas as large as those of the geographic oceanbasins. While the latter suggests a strong influence of theregions 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 ). Developed to investi-gate transition pathways in complex nonlinear stochasticsystems, TPT provides a statistical characterization ofthe ensemble of “reactive” trajectories, namely, pieces oftrajectories along which direct transitions between twosets A and B in phase space take place. The TPT termi-nology is borrowed from statistical mechanics and physi-cal chemistry, for which TPT was originally developed tostudy chemical reactions from reactants A to products B ,as an improvement for earlier approaches such as transi-tion state theory and transition path sampling. Sincethen, the TPT framework has also been applied to study- a r X i v : . [ phy s i c s . a o - ph ] S e p ing molecular conformation changes and transitionsin climate models. We here present, to the best ofour knowledge, the first oceanographic application.By constructing a Markov chain for debris motion andthen identifying coastline boxes in the ocean coveringwith reactant states A , and boxes in several ocean loca-tions including the subtropical gyres with product states B , we use TPT to infer pollution pathways in the globalocean. The Markov chain model accounts for the abilityof coastal boxes (states) to “pollute the oceans.” Thisinvolves adding an artificial state to the chain where alloutflow goes in and all inflow comes from (in an mannerthat differs from prior approaches ). By setting A toa single garbage patch and B as the union of the othergarbage patches, we can also assess the strength of thepatches.The rest of the paper is organized as follows. Theergodic-theory setup for closed systems is presented inSec. II. An adaptation of the theory for open systemsin discussed in Sec. III. The main results of TPT are re-viewed in Sec. IV, both for closed systems (Sec. IV A) andan 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 entailscoping with a number of issues, previously not encoun-tered, partially addressed, or overlooked, these in-clude: zonal connectivity; spurious communication be-tween ocean basins; and nonobserved communication; aswell as incorporating pollution sources near the coast.In Sec. VI time-asymptotic aspects of the chain dynam-ics are investigated, suggesting prospects for garbagepatches yet to be directly observed. The TPT analysisis applied in Sec. VII. This reveals pollution routes intothe garbage patches, which represent alternative targetsfor ocean cleanup efforts. Finally, a summary and theconclusions of the paper are presented in Sec. VIII.
II. SETUP FOR CLOSED DYNAMICAL SYSTEMS
Let us assume that floating debris trajectories are de-scribed by a time-homogeneous stochastic process in con-tinuous space
X ⊂ R and observed at discrete times nT , n ∈ Z . Its transition probabilities are controlled by astochastic kernel K ( x, y ) ≥ such that (cid:82) X K ( x, y ) y. = 1 for all x in phase space X , representing the world oceanbasin. The stochastic kernel is time-independent sincethe time-homogeneity of the process implies that therules 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 ) ≥ , (cid:82) X f ( x ) x. = 1 , describing the distribution of the randomposition X nT at any time nT evolves to the distribution P f ( y ) := (cid:90) X K ( x, y ) f ( x ) x. (1) at time ( n + 1) T , which defines a Markov operator P : L ( X ) (cid:9) generally known as a transfer operator . To infer the action of P on a discretized space onecan use a Galerkin projection referred to as Ulam’smethod. This consists of covering the phase space X with N connected boxes { B i } i ∈ S , S := { , . . . , N } ⊂ Z + , disjoint up to zero-measure intersections, and pro-jecting functions in L ( X ) onto the finite-dimensionalspace spanned by indicator functions on the boxes V N :=span (cid:8) Bi ( x )area( B i ) (cid:9) i ∈ S where A ( x ) = 1 if x ∈ A and 0 oth-erwise. The discrete action of P on V N is describedby a matrix P = ( P ij ) i,j ∈ S ∈ R N × N called a tran-sition matrix . The transition matrix results from theprojection P ij := Pr( X ( n +1) T ∈ B j | X nT ∈ B i )= 1area( B i ) (cid:90) B i (cid:90) B j K ( x, y ) x. y. (2)and describes the proportion of probability mass in B i that flows to B j during T . If one is provided with a largeset of observations x and x T of X and X T , respectively,then (2) can be estimated via counting the transitions inthe observed data, viz., P ij = C ij (cid:80) k ∈ S C ik , C ij := { x ∈ B i , x T ∈ B j } . (3)Note that (cid:80) j ∈ S P ij = 1 for all i ∈ S , so P is a row-stochastic matrix that defines a Markov chain on boxes,which represent the states of the chain.
The evolutionof the discrete representation of f ( x ) , i.e., the probabilityvector f = ( f i ) i ∈ S , (cid:80) i ∈ S f i = 1 , is calculated under left multiplication, i.e., f (cid:55)→ f P, (4)as it follows by noting that Pr( X ( n +1) T ∈ B j ) = (cid:80) i ∈ S Pr( X ( n +1) T ∈ B j , X nT ∈ B i ) = (cid:80) i ∈ S Pr( X nT ∈ B i ) P ij . In this paper, whenever we multiply vectors bymatrices, we assume that the vector takes the appropriateform of a row or column vector for the given operation.Because P is stochastic, = (1 , . . . , is a right eigen-vector with eigenvalue λ = 1 , i.e., P = . The eigen-value λ = 1 is the largest eigenvalue of P . The associatedpotentially nonunique left eigenvector p = ( p i ) i ∈ S is in-variant, because p P = p and can be chosen component-wise nonnegative (by the Perron–Frobenius theorem).We call P irreducible (or ergodic ) if for all i, j ∈ S thereexists n ij ∈ Z +0 \ {∞} such that ( P n ij ) ij > . To wit, allstates of an irreducible Markov chain communicate, theeigenvalue λ = 1 is simple, and the corresponding lefteigenvector p is strictly positive. We call P aperiodic (or mixing ) if there exists i ∈ S such that gcd { n ∈ Z +0 :( P n ) ii > } = 1 . No state of an aperiodic Markov chainis visited cyclically.If P is ergodic and mixing, then p , normalized to aprobability vector ( (cid:80) i ∈ S p i = 1 ), satisfies < p = p P =lim n ↑∞ f P n for any probability vector f . We call p aninvariant limiting probability vector or stationary distri-bution .We adopt the traditional notation with { X t } t ∈ Z in-stead of { X nT } n ∈ Z and write, for instance, P ij =Pr( X t +1 = j | X t = i ) , when this simplifies the nota-tion. In what follows we will assume that P is both er-godic and mixing, and the system is in stationarity, i.e., Pr( X t ∈ B i ) = p i for all t ∈ Z .The Markov chain model we will deduce from data inSec. V A is, however, open, thus not ergodic. For this rea-son, 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 andback into it. This can happen for instance when thedomain of interest is a subregion of the closed world oceandomain X or when trajectory data are only available in asubregion of X . Other possibilities include poor samplingof X , weak communication within, or the situation wedescribe in Sec. V. In every case the resulting dynamicalsystem 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 adiscrete-time-homogeneous Markov chain on a box cov-ering { B i } i ∈ O of the ocean domain X but the probabilityto transition from one box with index i ∈ O to anywhereelse in the domain O is no longer strictly since proba-bility mass can flow out of the domain. We denote thetransition matrix on the open domain by P O with en-tries given by P Oij := Pr( X t +1 = j | X t = i ) for i, j ∈ O .Since the rows of P O no longer have to add up to one, P O represents a substochastic matrix.We assume that a larger domain S ⊃ O exists on whichthe dynamics are closed, i.e., the transition matrix P onbox entries i, j ∈ S is stochastic. Furthermore, when wesay that the dynamics on the open domain is stationary,we actually mean that the dynamics on the larger, closeddomain is stationary with distribution p = ( p i ) i ∈ S , whilewe denote the restriction to the open domain by p | O =( p i ) i ∈ O .For further analysis it is often useful to artificially closethe open system. From the closure of P O , we can, forinstance, get an estimate of p | O . Closing P O can be doneby appending to O a state ω , which we will call two-waynirvana 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 notin O are lumped together, this restricted dynamics shouldbe consistent with the original one under the assumptionof well-mixedness between exit from O and reentry intoit. For simplicity of notation, we will denote the singleton { ω } also by ω and refer to it too as the two-way nirvanastate.The resulting transition matrix on O ∪ ω reads (possibly overloading the notation by denoting it by P again) P = (cid:18) P O P O → ω P ω → O (cid:19) (5)where P O → ω := (cid:0) − (cid:80) j ∈ O P Oij (cid:1) i ∈ O (understood as acolumn vector) gives the outflow from O to ω and P ω → O is a (row) vector that gives the inflow and has to be aprobability vector. Note that the matrix P is stochastic (cid:80) j ∈ O ∪ ω P ij = 1 for all i ∈ O ∪ ω and as such constitutesa closed dynamical system .When no information about the reentry is available,e.g., because data outside the open domain of interest arenot available, a possible choice for P ω → O is to redis-tribute according to the quasistationary distribution of P O . Lünsmann and Kantz alternatively use contouradvection for estimating the transition probabilities be-tween boxes. Without adding a nirvana state, Froyland,Pollett, and Stuart immediately redistribute the out-flow back into the system. Here we redistribute in sucha way that accounts for ocean pollution, as we describein Sec. V.In the next section we will see how to study transitionsbetween A and B (subsets of O ) in both the cases wherei) the domain is closed, i.e., O = S , and where ii) pathsonly traverse the open domain O (cid:40) S . In the lattercase, for the TPT computations only knowledge of P O and the estimate of the stationary density on the opencomputational domain p | O is necessary. IV. TRANSITION PATH THEORYA. TPT for closed systems
Motivated by a desire to understand rare events suchas transformations involved in chemical reactions, TPTprovides a rigorous approach to study transitions from aset A ⊂ S to another, disjoint set B ⊂ S of a Markovchain. The results presented below pertain to time-homogeneous (i.e., autonomous) chains; extensionsto the nonautonomous case have been recently derived, but they are beyond the scope of this paper. Tradition-ally, source set A is thought to be formed by reactant states, while target set B of product states. Thus transi-tions from A to B are referred to as reaction events , whilethe pieces of trajectories running from A to B without go-ing back to A or going through B in between are knownas reactive trajectories , which are the focus of TPT (Fig.1).The main tools of TPT are the forward and back-ward committor probabilities giving the probability of arandom walker to hit B before A , in either forward orbackward time. The committor probabilities are usedto express various statistics of the ensemble of reactivetrajectories: i) the density of reactive trajectories , whichprovides information about the bottlenecks during thetransitions; ii) the current of reactive trajectories indi-cating the most likely transition channels; iii) the rate of A BS
FIG. 1. Given a Markov chain taking values on S , the cartoonshows in red the reactive pieces of a trajectory connectingdisjoint sets A, B ⊂ S . reactive trajectories leaving A or entering B ; and iv) the mean duration of reactive trajectories. We will introducethese in the following. Recall that we assume the chainto be stationary with distribution p .The first entrance time of a set S ⊂ S is the stoppingtime random variable defined as τ + S := inf { t ≥ X t ∈ S} (6)where inf ∅ := ∞ . The forward committor q + := ( q + i ) i ∈ S gives the probability that a trajectory starting in i ∈ S first enters B , not A , i.e., q + i := Pr( τ + B < τ + A | X = i ) . (7)Note that q + i ∈ A = 0 while q + i ∈ B = 1 . For i ∈ C := S \ ( A ∪ B ) , one has that q + i = (cid:88) j ∈ S P ij q + j . (8)The solution to this algebraic system is unique due to theirreducibility of P , and in matrix notation expressed as q + | C = (cid:0) Id | C |×| C | − P | C (cid:1) − P | C,B | B |× , q + | A = | A |× , q + | B = | B |× , (9)where | S denotes the restriction on indices in S , while | S , S (cid:48) gives the restriction to rows corresponding to S andcolumns of S (cid:48) , if S = S (cid:48) we shorten this to | S .The last exit time , in turn, is defined by τ −S := sup { t ≤ X t ∈ S} (10)where sup ∅ := −∞ , which is a stopping time, but for the time-reversed chain { X − t } t ∈ Z that traverses the originalMarkov chain backwards in time, i.e. X − t := X − t . The reversed chain’s transition matrix, P − = ( P − ij ) i,j ∈ S isgiven by P − ij = Pr( X t = j | X t +1 = i ) = p j p i P ji , (11)since the chain is assumed to be in stationarity. Thetime-reversed transition matrix P − is ergodic and mix-ing, and has the same stationary distribution p as P .The backward committor q − := ( q − i ) i ∈ S gives the proba-bility that a trajectory starting in i ∈ S last exits A , not B : q − i := Pr( τ − A > τ − B | X = i ) . (12)In this case, q − i = (cid:88) j ∈ S P − ij q − j (13)for i ∈ C , subject to q − i ∈ B = 0 and q − i ∈ A = 1 . The(unique) solution in matrix notation, q − | C = (cid:0) Id | C |×| C | − P − | C (cid:1) − P − | C,A | A |× , q − | A = | A |× , q − | B = | B |× . (14)A particular situation arises in the special case when thechain is reversible, namely, when p i P ij = p j P ji or, equiv-alently, P − = P . In such a case, q − = − q + .The committors contain information that enable thecomputation of various transition statistics. The distri-bution of reactive trajectories µ AB = ( µ ABi ) i ∈ S , definedas the joint probability that the chain is in state i while transitioning from A to B , viz., µ ABi := Pr( X = i, τ − A > τ − B , τ + B < τ + A ) , (15)tells us where reactive trajectories spend most of theirtime. Note that µ ABi ∈ A ∪ B = 0 . The distribution of reac-tive trajectories is computable from the committor prob-abilities and the stationary distribution, µ ABi = q − i p i q + i . (16)A density of reactive trajectories ˆ µ AB = (ˆ µ ABi ) i ∈ S isobtained by normalizing µ ABi by the probability to bereactive Z AB := (cid:88) j ∈ C µ ABj = Pr( τ − A > τ − B , τ + B < τ + A ) , (17)as it follows from the law of total probability. The resultis ˆ µ ABi := µ ABi Z AB = Pr( X = i | τ − A > τ − B , τ + B < τ + A ) , (18)i.e., the probability of being in state i conditioned onbeing already on a reactive path from A to B .The current (or flux) of reactive trajectories f AB =( f ABij ) i,j ∈ S gives the average flux of trajectories goingthrough i and j at two consecutive times while on theirway from A to B : f ABij := Pr( X = i, X = j, τ − A > τ − B , τ + B < τ + A ) , (19)which is computable as f ABij = q − i p i P ij q + j . (20)Note that the reactive current can include direct tran-sitions from i ∈ A to j ∈ B , which are not accountedfor in the corresponding reactive distribution as it onlyconsiders transitions passing through C .To eliminate detours of reactive currents, one intro-duces the effective current of reactive trajectories f + =( f + ij ) i,j ∈ S , which gives the net amount of reactive currentgoing through i and j consecutively, viz., f + ij := max (cid:8) f ABij − f ABji , (cid:9) . (21)To visualize f + on a flow domain covered by boxes { B i } i ∈ S , one usually depicts the magnitude and the di-rection of the effective current out of each i , i.e., to each i one attaches the vector (cid:80) j (cid:54) = i f + ij e ij , where e ij is the unitvector pointing from the center of box B i to the centerof B j . There also exists a flow decomposition algorithmfor extracting the dominant transition paths from f + . The rate of transitions leaving A or departure rate isdefined as the probability per time step of a reactive tra-jectory to leave A , i.e., k A → := Pr( X ∈ A, τ + B < τ + A ) = (cid:88) i ∈ A,j ∈ S f ABij (22)and can be computed by summing up the reactive fluxthat exits A . In turn, the rate of transitions entering B or arrival rate is defined as the probability per time stepof a reactive trajectory to enter B : k B ← := Pr( X ∈ B, τ − A > τ − B ) = (cid:88) i ∈ S,j ∈ B f ABij . (23)By a simple calculation, it can be shown that summingthe reactive current out of A , (cid:80) i ∈ A,j ∈ S f ABij , is equal toaggregating the reactive current into B , (cid:80) i ∈ S,j ∈ B f ABij ,thus k A → = k B ← =: k AB . (24)To better interpret the transition rate k AB , we give twomeanings. Consider an infinite p -distributed ensembleof random walkers in our domain, then at any time theproportion of random walkers that are exiting A while ontheir way to B (or equivalently, entering B when cominglast from A ) is given by k AB . Now, on the other hand,consider only one random walker in the system, then k AB can be interpreted as a frequency, i.e., the random walkerexits A on average every ( k AB ) − -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 dis-connected set, it is insightful to further decompose thetransition rate k B ← = (cid:88) B n ⊂ B k B n ← (25)into the individual arrival rates into disjoint subsets B n that together give B = ∪ n B n : k B n ← = Pr( X ∈ B n , τ − A > τ − B ) = (cid:88) i ∈ S,j ∈ B n f ABij . (26)The same can also be done for decomposing k A → .Finally, dividing the probability of being reactive bythe discrete transition rate, t AB := Z AB k AB , (27)gives the expected duration of a transition from A to B . We close this section with a remark on comparing prob-abilistic computations with counting. Ergodicity of thechain implies that the objects in TPT can be approxi-mated by “counting” transition events of one sufficientlylong trajectory, and this approximation converges almostsurely as the length of the trajectory tends to infinity.
For instance, the forward committor q + i of any state i isapproximated by the fraction of all visits of the chain tostate i after which the chain directly transitioned to B without hitting A first. All other quantities consideredhere can be similarly approximated. As we intend toapply TPT to a chain extracted from drifter trajectorydata, one might wonder whether this level of sophistica-tion is necessary to our ends or whether one could simplydo an approximation by counting. The answer lies in thefeatures of the data. One would need sufficiently manydrifter trajectories that are sufficiently long to resolvethe transition statistics, and that are also spread accord-ing to the right distribution. None of these requirementsare met, and the best one can do is to “concatenate” thedrifter information into a Markov chain, as it will be donein Sec. V below. B. TPT for open domains
To apply TPT to open dynamical systems on O , amodification from the standard setting as reviewed inSec. IV A is needed. Adding the state ω to O closesthe system artificially (as in Sec. III) but we are stillonly interested in the transitions from A ⊂ O to B ⊂ O that stay in O during the transition. Thus the reactivetrajectories we consider go from A to B without passing A , B or ω during the transition. If we were to apply theusual TPT on the artificially closed system we would alsoobserve artificial transitions via the added state ω .In order to compute the statistics of the reactive tra-jectories from A to B only through O we look at slightlydifferent committors. Namely, the forward committornow gives the probability to next transition to B ratherthan to A or outside of O when starting in state i , i.e., q + i := Pr( τ + B < τ + A ∪ ω | X = i ) , (28)while the backward committor gives the probability tohave last come from A , not B ∪ ωq − i := Pr( τ − A > τ − B ∪ ω | X = i ) . (29)In that way the product of forward and backward com-mittors becomes the probability when initially in i tohave last come from A and next go to B while not pass-ing through A , B or ω in between.By definition, the forward committor is q + i = 0 for i ∈ A ∪ ω and for i ∈ B , while in the transition region C := O \ ( A ∪ B ) it satisfies q + i = (cid:88) j ∈ O ∪ ω P ij q + j = (cid:88) j ∈ O P ij q + j + P iω q + ω = (cid:88) j ∈ O P Oij q + j (30)since q + ω = 0 and P on entries of O reduces to P O .The backward committor q − i = 0 for i ∈ B and 1 for i ∈ A ∪ ω , while, by a similar reasoning as above, itsatisfies q − i = (cid:88) j ∈ O P O, − ij q − j (31)for i ∈ C , where P O, − is the restriction of the backward-in-time transition matrix P − to O and has entries P O, − ij = p j p i P Oji for i, j ∈ O .Therefore, system (9) remains the same with the re-placement of P with P O and A with A ∪ ω . In turn,system (14) remains the same with the replacement of P − with P O, − and B with B ∪ ω .The rest of the formulae in Sec. IV A are not changedexcept that the committors are now given as above. Animportant observation, however, is that µ ABi = 0 for i = ω and f ABij = 0 for i, j = ω . Thus only their valueson O are of interest, where P can be replaced by P O and p can be substituted by its restriction to O , p | O . Also,as the rate and mean transition time of reactive trajec-tories are derived from the density and current, they arecomputable solely from P O and p | O .This version of TPT for open dynamics, can also beuseful in other settings, e.g., when one wants to studytransitions 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 forthe dynamics of a single plastic debris piece that entersthe ocean at the coast with a probability reflecting ob-served levels of mismanaged plastic waste in near coastalcommunities. From the coast, the debris piece traverses the ocean, possibly passing and staying for long timesnear garbage patches. Its motion is fitted using satellite-tracked drifter trajectories; cf. Sec. V A. Whenever a de-bris piece beaches somewhere, we reinject it again nextto the coast. The coastal injection and beaching is de-scribed in Sec. V B below.In that way we will model the behavior of a genericplastic debris piece in the ocean by a stationary ergodicMarkov chain. Of course, there is a huge amount of plas-tic debris in the ocean, each day growing in number. Butwe are not interested in modeling the change in plasticconcentration in the ocean. Rather, our interest lies inunderstanding the routes of plastic waste from the coaststo the garbage patches by means of a TPT analysis. Thisdistinction is elaborated on in Sec. V C.
A. Preparation of P from drifter trajectory data As anticipated, to formulate the Markov chain for ma-rine debris motion we use drifter trajectory data, takenfrom the NOAA Global Drifter Program. Satellite-tracked by the
Argos system or GPS (Global PositioningSystem), the drifters from this database have a spher-ical surface float with a 15-m-long holey-sock drogueattached. They are engineered to resist wind slippageand wave-induced drift, and hence to follow water mo-tion as close as possible. We therefore only considertrajectory portions during which the drifter’s drogue hasbeen lost, which can be expected to provide a more fairrepresentation of floating marine debris motion. The basic procedure to construct the transition matrix P , defined in (3), is as follows. We first interpolate theavailable undrogued drifter trajectories daily and formtwo arrays, one representing positions at any instant oftime over 1992–2019 ( x ) and another one representingtheir images ( x T ) after T = 5 d. Here we are assumingthat the ocean motion did not change considerably overthe last 30 yr such that the transition matrix P fromthis data set is still a good representation of the “averageocean motion.”We then define the box covering { B i } i ∈ S by lying downon the world ocean domain a grid of roughly 3 ◦ width(due the planet’s curvature the area of the boxes is notfixed, varying from 100–10000 km , but this is inconse-quential in the definition of the vector space V N , nor-malized by box area). The entries of P are finally es-timated via counting according to (3). As in previouswork the transition time T is chosen long enoughto guarantee negligible memory into the past and suffi-cient communication among boxes, made large enoughto maximize sampling. The simple Markovianity test λ ( P ( nT )) = λ ( P ( T )) n is passed well up to n = 4 .There are additional aspects, not encountered, par-tially addressed or overlooked earlier, which must becoped with to make P meaningful.1. Zonal connectivity . This is addressed by iden-tifying and continuating trajectories crossing theantimeridian connecting the eastern and westernhemispheres.2.
Spurious communication between ocean basins .This situation occurs where ocean basins are sepa-rated by narrow land masses. The situations thatconcern us are the Panama Isthmus separating thePacific and Atlantic Ocean basins, and also theMaritime Continent separating the Pacific and In-dian Oceans. Neither the undrogued drifters con-sidered nor drogued drifters analyzed earlier re-veal connectivity between the Pacific and IndianOceans through the various straits and passages inthat region, which might seem at odds with thepresence of the Indonesian Throughflow, partic-ularly for the drifters drogued at 15 m. However,this takes place mainly within the thermocline layer(50–200 m), which is less correlated with localwind flow that quite strongly affects the undrogueddrifters and also the drogued drifters, albeit to alesser extent. To avoid spurious communication be-tween the basins we proceed as follows. Let B k bea box spanning portions of for instance the PacificOcean and Atlantic Ocean (Caribbean Sea). De-note B PO k and B AO k the portions of B k lying onthe Pacific Ocean and Atlantic Ocean sides, re-spectively. In computing transitions between B k and other boxes we only consider those from orinto B PO k or B AO k depending on which one makesthe largest number of transitions. This guaranteesthat P kj > and P jk > exclusively for j ∈ S such that B j is either in the Pacific Ocean or theAtlantic Ocean.3. Nonobserved communication . A prominent exam-ple of this is the communication between the At-lantic Ocean and Mediterranean Sea. Dependingon the size of the boxes B i , a connection might ex-ists through the Gibraltar Strait, even thought inreality no drifter is seen to traverse it (in any di-rection). We resolve this situation by excluding theMediterranean Sea domain from consideration.4. Weak communication . We enable as much com-munication as possible along the chain by restrict-ing the chain to the largest strongly communicatingclass of states. This is done by applying the Tarjanalgorithm on the directed graph equivalent to theMarkov chain. This procedure excludes boxes fromthe partition. Among those boxes are 22 poorlysampled coastal boxes, mainly in the Kara Sea ofthe Arctic Ocean and the Seas of Indonesia, withtrajectories flowing in, but not flowing out in thenext step. Let O be the ordered set of box indices inthe largest class of strongly connected boxes. Usingthe notation in Sec. III, we call P O the substochas-tic transition matrix characterizing this open sys-tem. The Markov chain is now substochastic, sinceby the exclusion of boxes it is no longer ensuredthat probability mass is conserved. B. Pollution-aware model derivation
To formulate our Markov-chain model for ocean pol-lution, we leverage the possibility that marine debris getstuck at shorelines. This creates additional outflow of thesystem that must be compensated for, which we chooseto do in such a way as to model ocean pollution at thecoasts.Specifically, let (cid:96) : O → [0 , be a land fraction func-tion giving the ratio between land area and total boxarea. Namely, < (cid:96) ( i ) < for i ∈ L ⊂ O correspondingto boxes filled with some portion of land (or ice) (Fig. 2,top panel) and (cid:96) ( i ) = 0 otherwise. We then follow Miron et al. and replace P Oij ← (cid:0) − α(cid:96) ( i ) (cid:1) P Oij , ≤ α ≤ (32)for all i, j ∈ O . To wit, only a fraction of the probabilitymass, proportional to the amount of land covering box B i , is allowed to flow from i to j , the remaining probabil-ity mass is assumed to beach and flows out of the system.The factor α , not considered in Miron et al. , was in-cluded to enable consistency with observations. Whilewe have performed optimizations of no kind, we havefound that α = produces results most consistent withthem. If α = 1 (as in Miron et al. ) then the so-calledGreat Pacific Garbage Patch in the North Pacific sub-tropical gyre is not revealed as intense as observationsindicate. However, transition channels into this patchand patches in the other subtropical gyres are not sensi-tive to the specific α -value assumed, as we show in thesupplementary material.To deal with the created substochasticity by a closureof the system, we augment the chain by one artificialstate ω as in (5). All the outflow of the open system goesinto ω and we reinject the probability mass from ω to O through coastal boxes according to plastic waste inputfrom land into the ocean, viz., P ω → Oi = W i (cid:80) i ∈ L W i if i ∈ L, otherwise . (33)Here W i is the mass of mismanaged plastic waste in B i , i ∈ L , as inferred from estimates made in 2010 forpopulations living within 50 km of the coastline. Thisis shown in percentage of the total mass in the bottompanel of Fig. 2; note that only inhabited coastal boxesfor which estimates are available are shown. We denotethe transition matrix on the closed domain by P , butit should not be confused with the transition matrix P from the above Sec. V A which has a different domain andentries. The two-way nirvana state ω compensates for thesubstochasticity of P O = P | O by sending back into thechain any imbalances through the land states distributedaccording to the ability of such states to “pollute theoceans” as inferred by their share of the global land-basedplastic waste entering the ocean through them. It mustbe realized that in this statistical model debris mass is FIG. 2. (top panel) Fraction of land (or ice) filling coastalboxes of the surface world ocean partition (black). (bottompanel) Percentage of share of global mismanaged littered orinadequately disposed plastic waste estimated in 2010 for pop-ulations living within 50 km of the coastline. neither created nor destroyed. In other words, the modelassumes that the world ocean is polluted by plastic ata certain level, and that the ocean currents and windsredistribute the existing pollutants within ocean basins.If beaching occurs, then the pollutants are returned backto the ocean in an equal quantity simulating mismanagedplastic waste loading from land runoff.
C. Physical interpretation of the model
We model the distribution of garbage input pertime unit to the oceans by a time-independent vector P ω → O =: r ∈ R | O | . Each entry of r accounts for theprobability per time unit of injecting a garbage particleinto the corresponding box. Thus r is supported on thecoastal (land) boxes, i.e., r i = 0 for i / ∈ L , and r is aprobability vector, i.e., (cid:80) i r i = 1 .Then the total accumulated garbage mass distribution in the oceans is time-asymptotically going to be ∞ (cid:88) k =0 r ( P O ) k = r (Id − P O ) − . (34)Recall that P O is assumed to be irreducible, thus Id − P O is invertible. Equation (34) gives the mass distributionof debris particles entered over an infinite time frame,thus it does not need to be a probability vector. It is thelimiting (saturated) mass distribution of pollution mea-sured in the units dictated by r . If we would like to know the relative distribution of garbage that has accumulatedover time, we would norm this vector to a probabilityvector.Now, it turns out that the very same long-term distri-bution is modeled by our “recirculating” Markov chain.With a := P O → ω being the vector of absorption prob-abilities from the boxes into nirvana (the outflow), thestationary distribution of our chain satisfies (cid:0) p | O ρ (cid:1) (cid:18) P O ar (cid:19) = (cid:0) p | O ρ (cid:1) , (35)with stationary vector p | O on indices of O and scalar ρ giving the stationary weight of ω . The set of equationscorresponding to the boxes in O read as p | O P O + ρ r = p | O , or, after rearrangement, p | O (Id − P O ) = ρ r . (36)Since ρ is scalar, this readily means that p | O ∝ r (Id − P O ) − , which equals the asymptotic mass distri-bution (34) from above. In summary, our stationaryMarkov chain constructed with reinjection is the statisti-cal equivalent of garbage motion in the ocean, based onthe limiting garbage distribution. VI. LONG-TIME ASYMPTOTICS
By design, the proposed transition matrix P for ma-rine debris pollution has a single maximal communicat-ing class of the states, implying irreducibility for P andergodicity for the dynamics. Furthermore, direct push-forward (i.e., evolution under left multiplication by P ) ofan arbitrary probability vector reveals convergence to thedominant left eigenvector p (the chain is also aperiodic),which is invariant and also limiting, and hence repre-sents a stationary distribution. The top panel of Fig. 3shows p > restricted to O , viz., the set of boxes ofthe world ocean partition where the dynamics are open.The middle and bottom panels show, restricted to O ,the distribution after 1 and 10 yr of evolution under leftmultiplication by P of ω , respectively. Note that p | O lo-cally maximizes in the subtropical gyres, quite evidentlyin the eastern side of the North Pacific gyre. In mostof the Indian Ocean p | O reveals several well-spread localmaxima consistent with a predominantly uniform distri-bution. The exception is the Bay of Bengal, where p | O shows more clear sings of local maximization. An addi-tional local maximum of p | O is seen in the Gulf of Guineasouth of West Africa. The several local maxima of p | O identified are indicated by the red boxes in the top panelof Fig. 3. The Indian Ocean location corresponds to itslocal maximum inside the subtropical gyre.The structure of p | O suggests garbage patches inthe subtropical gyres of the Atlantic and PacificOceans consistent with in-situ microplastic concentra-tion observations. Previous analyses of drifter datarevealed these patches too, albeit from direct evolution
FIG. 3. (top panel) Restricted to the set O of boxes coveringthe physical ocean domain where the dynamics are open, thestationary distribution p of the closed dynamics representedby the transition matrix P for marine debris pollution. Notethat p locally maximizes inside the subtropical gyres, which,at the same time happen to develop great patches. Indicatedby the red boxes are these (and additional; cf. text for details)local maxima of p | O . (middle panel) Restricted to O , distri-bution after 1-yr evolution under left multiplication by P of aprobability density (vector) with support on the virtual nir-vana state included to close the system. (bottom panel) Asin the middle panel, but after 10 yr. of probability densities. In particular, van Sebille, Eng-land, and Froyland argued that the North Pacific patchshould be the main attractor of global marine debris,in agreement with direct observational evidence ofthe Great Pacific Garbage Patch. Our pollution-awaremodel produces consistent results. This can be antici-pated from p | O acquiring larger values in the North Pa-cific gyre than in the other subtropical gyres, and alsofrom direct pushforward of ω and subsequent restric-tion of the evolved density to the boxes where p | O lo-cally maximizes in the subtropical gyres (Fig. 4). (Itshould be noted too that the structure of p | O in the FIG. 4. Evolution of ω under left multiplication by P re-stricted to the boxes where p | O locally maximizes in the sub-tropical gyres. North Pacific suggests a garbage patch, albeit weaker,in the western side of the basin in agreement with fieldsampling. ) The relative weakness of the Indian Oceangarbage patch attributed to unique oceanic and atmo-spheric dynamics in the region is consistent with theresults from our Markov-chain model for ocean pollu-tion too. There the stationary distribution p | O does notreveal a clear local maximum (Fig. 3), and the directpushforward of ω identifies the Indian Ocean gyre asthe less attracting of all the subtropical gyres (Fig. 4).Exactly where the garbage patches are located is deter-mined by wind-induced Ekman and wave-induced Stokesdrift effects mediated by the inertia (i.e., buoyancy andsize) of the floating debris pieces. Indeed, thenumerical simulations of inertial particles by Beron-Vera,Olascoaga, and Lumpkin do not reveal signs of accumu-lation in the Indian Ocean gyre as clear as in the othergyres.In van Sebille, England, and Froyland the authorssuggest the possibility of a rather strong garbage patchin the Barents Sea in the Arctic Ocean, possibly con-strained by slow surface convergence due to deep-waterformation. While the authors noted that this patch mightbe an artifact of drifters becoming grounded in the (sea-sonal) sea-ice, observational support of plastic accumu-lating in the region is emerging. However, the observedaccumulation represents a very small fraction (3%) of theglobal standing stock. Our pollution-aware model doesnot reveal a patch there, more consistent with this ob-servation.Our model suggests the occurrence of a patch in theGulf of Guinea, which seems to be supported only on nu-merical simulations. However, the the Gulf of Guineais identified as a mesopelagic niche with genomic char-acteristics than different than its surroundings. Thispatch remained elusive to earlier studies.
A likely ex-planation is the involvement in those earlier studies ofboth undrogued and drogued drifters, which unlike float-ing debris, are much less affected by inertial effects. However, a more recent study involving exclusively un-0drogued drifters did not reveal accumulation in the Gulfof Guinea time asymptotically.The structure of p | O also reveals that the Bay of Ben-gal has potential for holding a garbage patch. Highplastic concentration in the Bay of Bengal has been re-ported and attributed to loading from nearby land-basedsources. The occurrence of a garbage patch in the Bayof Bengal was also suggested recently from the analysisof undrogued drifter trajectory data. The pertinent question is how the garbage patches arefilled. We address this using TPT.
VII. REACTIVE DEBRIS PATHS
With the above in mind, we proceed to apply TPT tothe dynamics on the physical world ocean domain, wherereactive debris currents are sought to be unveiled. Theusual TPT (Sec. IV A) allows us to compute statisticsof the ensemble of reactive paths of marine debris intogarbage patches, with the help of TPT for open domains(Sec. IV B) we can study reactive paths between garbagepatches.
A. Pollution paths into garbage patches
To infer the pollution paths into garbage patches, wechoose the nirvana state ω as the source state A ofgarbage, and we identify the set of target states B withthe union of indices of boxes covering garbage patches asinferred by the regions where p | O tends to locally maxi-mize, which we have isolated above (cf. Fig. 3, top panel).We will denote G the garbage patch set. Although thedebris is reentering the ocean through the land boxes L ,choosing the source as A = ω is more reasonable, as itallows reactive debris trajectories to enter boxes in L andto flow on towards B . With the choice A = L we wouldhave excluded this possibility, which would have causeda notable impact on TPT computations, given the size ofour boxes. The effective currents of reactive trajectoriesresulting 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 coastalboxes L .We first note that the extent of the reactive currentsrunning into the subtropical gyre patches is in generallarger than those running into the near coastal patches.These indicates that the near coastal patches are mainlyfed from nearby land-based plastic waste sources. Anexception is the Bay of Bengal patch, which appears toaccumulate garbage from remote sources in the coastsof the Arabian Sea and even more remote ones in thecoasts of Indonesia. Particularly constrained seems to bepatch in the Gulf of Guinea, which is inferred to be filledwith plastic debris releases at the southern coasts of WestAfrica. A refined assessment of these mostly qualitativeconclusions is presented below. Continuing with the visual inspection of Fig. 5, forthe North Pacific patch TPT analysis infers a robustzonal eastward reactive channel into it straight out fromthe coasts of China. A good deal of the transporteddebris is inferred to travel back to the western side ofthe North Pacific basin, where the stationary distribu-tion p | O also tends to maximize. A pollution source forthe South Pacific patch is not restricted to the east coastof South America. Indeed, a westerly transition chan-nel originating in the Indian Ocean and the coasts ofNew Zealand is also identified. Two clear reactive pathsinto the North Atlantic are identified, one mainly com-ing from the southeastern coast of the United States andanother one coming from the northern coasts of WestAfrica. In turn, the South Atlantic patch is fed from de-bris transport from the Brazil–Malvinas Confluence andthe southern tip of Africa. A main carrier of pollution forthe Indian Ocean patch is the Agulhas Return Current.However, this pollution channel bifurcates a bit east ofthe patch’s longitude, where a branch originates to ulti-mately feed the South Pacific patch. Indeed, the patternof the currents near the Indian Ocean patch does not sug-gest as clear channels into it as into the patches in theother subtropical gyres. This seems consistent with thereported weakness of the Indian Ocean patch.It is important to note that while ocean currentsplay a dominant role in transporting debris, the reac-tive paths inferred by TPT do not resemble entirely themean surface-ocean currents. However, this is not unex-pected given the various mechanisms, noted above, con-trolling the motion of floating material beyond advec-tion by ocean currents. We stress again that TPT, byconstruction, highlights currents composed of only tra-jectories that go from A (source) to B (target), therebyexcluding information about currents that go from B to A , A to A , and B to B .The expected transition duration (27) is estimated tobe 2.6 yr from the coasts into the subtropical gyre patchesand the gulf and coastal sea patches. Note that this isthe mean time a reactive trajectory takes from being in-jected into the oceans to hit any of these patches. If weset B = Y , where Y ⊂ G is the set of indices corre-sponding to the subtropical gyre patches, then the meanduration is 5.6 yr, cf. (26) and (27). The expected dura-tions of individual transition paths into the North Pacific,South Pacific, North Atlantic, South Atlantic, and IndianOcean patches are 7.3, 8.6, 4.3, 4.0, and 4.2 yr, respec-tively. The mean durations of those into the patches inthe Bay of Bengal and the Gulf of Guinea are 0.6 and 0.2yr, respectively. The proximity to the coasts explain theshort mean durations of the latter transition channels.As for the transition channels into the subtropical gyrepatches, those into the South Pacific and North Atlanticpatches stand out as the overall slowest and fastest inthe class, respectively. These times to individual patchesrepresent the mean duration of those reactive trajectoriesthat first hit the set B = g ∈ G through the respectivepatch, i.e., transitions are direct and not through other1 FIG. 5. Inferred reactive probability currents of marine debris into garbage patches (red boxes). Black boxes indicate coastalboxes from which those currents emerge. patches, G \ g (which can be avoided by using TPT foropen dynamics).Additional insight is provided by the normalized dis-tribution of reactive trajectories ˆ µ AB , plotted in Fig. 6,showing where reactive trajectories spend most of thetime while on their way from source A to target B . Notethat the reactive trajectories tend to bottleneck over largeregions around the subtropical gyre patches except theIndian Ocean gyre patch. These regions measure the sizeof the patches. The bottleneck is particularly pretty in-tense in the North Pacific gyre. The reactive flows inthe Indian Ocean patch are not seen to spend as muchtime near the patch as near the other subtropical gyrepatches. This is consistent with it being a weak garbagepatch. Additional intense bottlenecks are observed toconcentrate in the Bay of Bengal.Further insight into the domain of influence of each in-dividual garbage patch g ∈ G , and thus into the locationson the coast where debris flows into them originate from,is offered by associating to each state i ∈ O the mostlikely patch g (target) to hit according to the probabilityin i to forward-commit to g , viz., q + i ( g ) = Pr( τ + g < τ + ω | X = i ) . (37)This way every box of the partition gets assigned to apatch, forming what we call a forward-committor-baseddynamical geography , which is shown in Fig. 7. Note thelarge influence exerted by the subtropical patches on theglobal transport of marine debris, particularly those inthe subtropical gyres whose provinces span the largestareas. Similar influence of the subtropical patches was FIG. 6. Probability density of reactive debris paths that indi-cating where debris bottlenecks in their way into the garbagepatches. inferred from spectral analysis applied on simulatedtrajectories and from direct evolutions using drifter tra-jectory data. The relatively large influence of the Bay ofBengal patch inferred from the visual inspection of thereactive currents into it is well framed by the geography.The provinces of the geography in Fig. 7 are coloredaccording to the mean residence time , defined as follows.Let Q ⊂ O be the box indices of a given province. Themean time it takes a trajectory initialized in i ∈ Q tomove out of Q and thus hit the complement of Q , h Qi := E ( τ + O ∪ ω \ Q | X = i ) , is given by the solution of the linearequation (cid:0) Id | Q |×| Q | − P | Q (cid:1) h Q = | Q |× , (38)2 FIG. 7. Forward-committor-based dynamical geography re-vealing domains of influence for the garbage patches with theprovinces colored according to residence time. where h Q = ( h Qi ) i ∈ Q . By taking the average of h Q withrespect to the stationary density p | Q we get the residencetime in Q , i.e., H Q := E ( τ + O ∪ ω \ Q | X ∈ Q ) = h Q · p | Q p | Q · | Q |× . (39)The longest residence time is 14.6 yrs, computed for theSouth Pacific province, whereas the shortest residencetimes are 0.7 and 0.3 yrs for the Bay of Bengal and Gulfof Guinea regions, respectively. The North Pacific Ocean,North Atlantic Ocean and South Atlantic Ocean subtrop-ical subtropical garbage patches all have comparable res-idence times that range between 7–7.5 yrs while the In-dian Ocean garbage patch has a much shorter residencetime of 1.8 years. B. Pollution paths out of subtropical garbage patches
The interconnectivity of the subtropical garbagepatches with respect to the amount of debris particlesthat are exchanged between patches is presented in Fig.8. More precisely, we compute the reactive flux from A = y ∈ Y ⊂ G to B = ( Y \ y ) ∪ ω , where Y is theset of subtropical gyre patches. Then, the proportions oftotal debris mass present in the ocean that flow per timestep out of A and make their way towards B are k AB =1 . × − , . × − , . × − , . × − , and . × − for A chosen as the North Pacific, South Pacific, NorthAtlantic, South Atlantic, and Indian Ocean patches, re-spectively. We can further decompose the transition ratefrom A to B into the sum of arrival rates into each indi-vidual patch b in B , k AB = k B ← = (cid:80) b ∈ B k b ← as in (26).For a fixed A , the arrival rates into each b ∈ B are shownin the rows of Fig. 8. Consistently, the “emission” from agarbage patch y recirculates almost completely through ω before reaching any other patch Y \ y , hence the muchhigher rates in the column corresponding to the nirvanastate. In addition, relatively high reactive rates betweenthe subtropical garbage patches of the southern hemi-sphere highlight an interconnection between the IndianOcean, the South Atlantic and the South Pacific patches. Specifically, the Southern Atlantic debris transit at highrate to the South Pacific and Indian Oceans and, sim-ilarly, the Indian Ocean debris transit at high rate tothe South Pacific and Atlantic Ocean garbage patches.Finally, the reactive rates from the North Atlantic gyreto any other subtropical garbage patches are negligible,confirming again that it has very little connection withother patches and debris that manage to escape it mostlikely end up on land or in ice. FIG. 8. Reactive rates from each subtropical gyre patch, pre-sented by row, into all other subtropical gyre patches and thenirvana state. The last row shows the rates from the nirvanastate ω into the subtropical gyre patches from the results pre-sented on Fig. 5. The last row of Fig. 8 shows transition rates from(26) corresponding to the currents into each subtropi-cal patches presented Fig. 5 with A chosen as the nir-vana state ω . As expected, the transition rates from ω to the subtropical garbage patches are orders of magni-tude higher than the transition rates between patches.Bearing in mind that those transition rates are very low,meaning that the transitions are unlikely, associated re-active currents are depicted in Figs. 9 and 10. These rep-resent potential pathways that marine debris might takeout of the gyres, for instance, in the event of unusuallystrong winds.Figure 9 presents the reactive currents from the sub-tropical gyre patches to the nirvana state ω , which corre-spond to the last column of Fig. 8. That is, we set A = Y (black squares) and B = ω (red squares are coastal bins i ∈ L where P O → ω > ). In general, debris out of thenorthern hemisphere patches have a larger probability ofbeaching than the southern hemisphere patches. In par-ticular, the reactive currents in the Indian Ocean followthe general path of debris from the search area of the in-famous Malaysia Airlines flight MH370 to the locations of3recovered debris on the coasts of Mauritius, Madagascar,Mozambique, Tanzania, and South Africa. In turn, Fig. 10 presents the reactive currents from asubtropical gyre patch to the union of all other subtrop-ical gyre patches. To place the focus on debris trajec-tories that stay in the ocean, we do not allow reactivepassages via ω . Thus we use TPT for open domains bysetting A = y (black square) and B = Y \ y (red squares)in (28) and (29). The reactive currents out of the In-dian Ocean patch are quite strong, in agreement withreports on its weak character. However, these are some-what weaker than those out of the South Atlantic patch.Note that both the Indian Ocean patch and the SouthAtlantic patch exchange debris with the South PacificOcean patch, as shown in Fig. 8, through the AntarcticCircumpolar Current. The currents that flow out of theNorth Pacific patch are much weaker, yet not as weakas those coming out of the North Atlantic patch. Thestrength of the currents out of the South Pacific patchranges in between the above.To quantify the above qualitative conclusions from theinspection of the transition channels, we computed thereactive rates from each y ∈ Y to Y \ y , telling us theamount of debris probability mass that flows out of y ∈ Y per time step and is on its direct way to Y \ y , equal to therow sums of Fig. 8 excluding the portion that goes intonirvana. The reactive rate (24) gives . × − , . × − , . × − , . × − , and . × − for the NorthPacific, South Pacific, North Atlantic, South Atlantic,and Indian Ocean patches, respectively, which confirmour qualitative assessments above. We note that eachof these rates is at least one order of magnitude smallerthan those reported at the very beginning of this section,except that of the South Atlantic, where it is merely afactor 5 weaker. This indicates that debris leaving theSouth Atlantic is most frequently finding its way to otherpatches.It must be noted that the above reactive rates do notsay anything about retention. They tell us which patch“emits” the most frequently such debris that finds its wayto another patch. A low rate does not need to mean thatdebris leaving A comes back to A since the debris canhit ω too before hitting B , as shown by the much higherreactive rates to the nirvana state in Fig. 8. In otherwords, a low rate should not be taken to mean the sameas high attraction. Thus, the reactive rate computationresults just described do not contradict those from directdensity evolution in Fig. 4, which had identified the NorthPacific patch as the most attracting of all. VIII. SUMMARY AND CONCLUSIONS
We have presented a novel application of transitionpath theory (TPT), here extended to open autonomousdynamical systems. The problem chosen was that ofpollution routes from possible coastline sources (reactivestates) into garbage patches in the global surface ocean
FIG. 9. Reactive currents from each subtropical gyre patchto the nirvana state ω . The source set is indicated in black ineach panel; the target sets are indicated in red. (product states).Undrogued drifter trajectories from NOAA GlobalDrifter Program were used to derive a Markov chainon which TPT was applied, as a model for the time-asymptotic dynamics of marine debris pollution. Mod-eling the probability of trajectories to beach as a func-tion of the fraction of land filling each coastal box of thecovering of the world ocean domain resulted in an opensystem, which was closed by sending the probability im-balance back into the chain according to the capacityof coastal boxes to “pollute the oceans” as measured byits share of global mismanaged plastic waste. Assuminga constant pollution rate, our time-homogeneous modelwas shown to be the statistical equivalent of a “saturated”(stationary) pollution redistribution dynamics.A high probability transition channel was identifiedconnecting the Great Pacific Garbage Patch with thecoasts of Eastern Asia, suggesting an important sourceof plastic pollution there. The weakness of the IndianOcean gyre as a trap of plastic debris was found con-sistent with transition paths not converging in the gyre.While the North Pacific subtropical gyre was found to bemost attracting consistent with earlier assessments, theSouth Pacific gyre stood out as the most enduring in thesense that the total reactive rate out of that gyre intoother gyres and the nirvana state resulted the smallestof all. The weakest of all the gyres in terms of its ca-pacity to trap and hold within plastic waste resulted tobe South Atlantic gyre. The gyres were found in generalweakly communicated. Indeed, in the event of anoma-lously intense winds a subtropical gyre is more likely toexport garbage out toward the coastlines than into an-other gyre.Our results, including prospects for garbage patchesyet to be directly and/or robustly observed, namely, theGulf of Guinea and the Bay of Bengal, have implica-tions for activities such as ocean cleanup as the revealedreactive pollution routes provide targets, alternative tothe great garbage patches themselves, to aim those ef-forts. Additional ocean applications of TPT are un-derway (e.g., using submerged float data and targeting4meridional overturning routes) and will be reported else-where. SUPPLEMENTARY MATERIAL
The supplementary material contains versions of Figs.4 and 5 assuming α = (Figs. S1 and S4, respectively), (Figs. S2 and S5), and (Figs. S3 and S6) in (32). ACKNOWLEDGMENTS R. Lumpkin and M. Pazos, “Measuring surface currents with Sur-face Velocity Program drifters: the instrument, its data and somerecent results,” in
Lagrangian Analysis and Prediction of Coastaland Ocean Dynamics , edited by A. Griffa, A. D. Kirwan, A. Mar-iano, T. Özgökmen, and T. Rossby (Cambridge University Press,2007) Chap. 2, pp. 39–67. E. van Sebille, E. H. England, and G. Froyland, “Origin, dynam-ics and evolution of ocean garbage patches from observed surfacedrifters,” Environ. Res. Lett. , 044040 (2012). A. N. Maximenko, J. Hafner, and P. Niiler, “Pathways of marinedebris derived from trajectories of Lagrangian drifters,” Mar. Pol-lut. Bull. , 51–62 (2012). 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, andC. M. Duarte, “Plastic debris in the open ocean,” Proc. Nat.Acad. Sci. USA , 10239–10244 (2014). F. J. Beron-Vera, M. J. Olascoaga, and R. Lumpkin, “Inertia-induced accumulation of flotsam in the subtropical gyres,” Geo-phys. Res. Lett. , 12228–12233 (2016). F. J. Beron-Vera, M. J. Olascoaga, and P. Miron, “Building aMaxey–Riley framework for surface ocean inertial particle dy-namics,” Phys. Fluids , 096602 (2019). F. J. Beron-Vera, “Nonlinear dynamics of inertial particles in theocean: From drifters and floats to marine debris and
Sargassum ,”Nonlin. Dyn. submitted, arXiv:2007.15638 (2020). A. Lasota and M. C. Mackey,
Chaos, Fractals and Noise:Stochastic Aspects of Dynamics , 2nd ed., Applied Mathemati-cal Sciences, Vol. 97 (Springer, New York, 1994). P. Brémaud,
Markov chains , Gibbs Fields Monte Carlo Simula-tion Queues, Texts in Applied Mathematics, Vol. 31 (Springer,New York, 1999). J. Norris,
Markov Chains (Cambridge University Press, 1998). M. Dellnitz and A. Hohmann, “A subdivision algorithm for thecomputation of unstable manifolds and global attractors,” Nu-merische Mathematik , 293–317 (1997). M. Dellnitz and O. Junge, “On the approximation of complicateddynamical behavior,” SIAM J. Numer. Anal. , 491–515 (1999). G. Froyland and M. Dellnitz, “Detecting and locating near-optimal almost-invariant sets and cycles,” SIAM J. Sci. Comput. , 1839–1863 (2003). P. Koltai,
Efficient approximation methods for the global long-term behavior of dynamical systems – Theory, algorithms andexamples , Ph.D. thesis, Technical University of Munich (2010). G. Froyland, R. M. Stuart, and E. van Sebille, “How well-connected is the surface of the global ocean?” Chaos , 033126(2014). E. Vanden-Eijnden, “Transition Path Theory,” Lect. Notes Phys. , 439–478 (2006). P. Metzner, C. Schütte, and E. Vanden-Eijnden, “Illustration oftransition path theory on a collection of simple examples,” J.Chem. Phys. , 084110 (2006). E. Weinan and E. Vanden-Eijnden, “Towards a theory of transi-tion paths,” J. Stat. Phys. , 503–623 (2006). P. Metzner, C. Schütte, and E. Vanden-Eijnden, “Transition paththeory for markov jump processes,” Multiscale Modeling & Sim-ulation , 1192–1219 (2009). E. Weinan and E. Vanden-Eijnden, “Transition-path theory andpath-finding algorithms for the study of rare events,” Annu. Rev.Phys. Chem. , 391–420 (2010). E. Wigner, “The transition state method,” Trans. Faraday Soc. , 29–41 (1938). L. Pratt, “A statistical method for identifying transition statesin high dimensional problems,” J. Chem. Phys. , 5045–5048(1986). F. Noé, C. Schütte, E. Vanden-Eijnden, L. Reich, and T. R.Weikl, “Constructing the equilibrium ensemble of folding path-ways from short off-equilibrium simulations,” Proceedings of theNational Academy of Sciences , 19011–19016 (2009). V. A. Voelz, G. R. Bowman, K. Beauchamp, and V. S. Pande,“Molecular simulation of ab initio protein folding for a millisecondfolder ntl9 (1- 39),” Journal of the American Chemical Society , 1526–1528 (2010). D. Lucente, S. Duffner, C. Herbert, J. Rolland, and F. Bouchet,“Machine learning of committor functions for predicting high im-pact climate events,” arXiv preprint arXiv:1910.11736 (2019). J. Finkel, D. S. Abbot, and J. Weare, “Path properties of atmo-spheric transitions: illustration with a low-order sudden strato-spheric warming model,” Journal of the Atmospheric Sciences , 2327–2347 (2020). G. Froyland, P. K. Pollett, and R. M. Stuart, “A closing schemefor finding almost-invariant sets in open dynamical systems,”Journal of Computational Dynamics , 135 (2014). B. Lünsmann and H. Kantz, “An extended transfer operator ap-proach to identify separatrices in open flows,” Chaos: An Inter-disciplinary Journal of Nonlinear Science , 053101 (2018). P. Miron, F. J. Beron-Vera, M. J. Olascoaga, J. Sheinbaum,P. Pérez-Brunius, and G. Froyland, “Lagrangian dynamical geog-raphy of the Gulf of Mexico,” Scientific Reports , 7021 (2017). P. Miron, F. J. Beron-Vera, M. J. Olascoaga, and P. Koltai,“Markov-chain-inspired search for MH370,” Chaos: An Interdis-ciplinary Journal of Nonlinear Science , 041105 (2019). P. Miron, F. J. Beron-Vera, M. J. Olascoaga, G. Froyland,P. Pérez-Brunius, and J. Sheinbaum, “Lagrangian geography ofthe deep Gulf of Mexico,” J. Phys. Oceanogr. , 269–290 (2019). M. J. Olascoaga, P. Miron, C. Paris, P. Pérez-Brunius, R. Pérez-Portela, R. H. Smith, and A. Vaz, “Connectivity of Pulley Ridgewith remote locations as inferred from satellite-tracked driftertrajectories,” Journal of Geophysical Research , 5742–5750(2018). F. J. Beron-Vera, N. Bodnariuk, M. Saraceno, M. J. Olascoaga,and C. Simionato, “Stability of the Malvinas Current,” Chaos ,013152 (2020). E. Morrison, A. Shipman, S. Shrestha, E. Squier, andK. Stack Whitney, “Evaluating The Ocean Cleanup, a marinedebris removal project in the North Pacific Gyre, using SWOTanalysis,” Case Studies in the Environment , 1–6 (2019). S. M. Ulam,
A Collection of Mathematical Problems , Intersciencetracts in pure and applied mathematics (Interscience, 1960). Z. Kovács and T. Tél, “Scaling in multifractals: Discretization ofan eigenvalue problem,” Phys. Rev. A , 4641–4646 (1989). L. Helfmann, E. R. Borrell, C. Schütte, and P. Kotai, “Ex-tending transition path theory: Periodically driven and finite-time dynamics,” J. Nonlinear Sci. doi.org/10.1007/s00332-020-09652-7 (2020). A. L. Sybrandy and P. P. Niiler, “WOCE/TOGA Lagrangiandrifter contruction manual,” Tech. Rep. SIO Reference 91/6(Scripps Institution of Oceanography, La Jolla, California, 1991). P. P. Niiler and J. D. Paduan, “Wind-driven Motions in the north-eastern Pacific as measured by Lagrangian drifters,” J. Phys.Oceanogr. , 2819–2830 (1995). R. Lumpkin, S. A. Grodsky, L. Centurioni, M.-H. Rio, J. A. Car-ton, and D. Lee, “Removing spurious low-frequency variability indrifter velocities,” J. Atm. Oce. Tech. , 353–360 (2012). M. J. Olascoaga, F. J. Beron-Vera, P. Miron, J. Triñanes, N. F.Putman, R. Lumpkin, and G. J. Goni, “Observation and quan-tification of inertial effects on the drift of floating objects at theocean surface,” Phys. Fluids , 026601 (2020). P. Miron, S. Medina, M. J. Olascaoaga, and F. J. Beron-Vera,“Laboratory verification of a Maxey–Riley theory for inertialocean dynamics,” Phys. Fluids , 071703 (2020). R. McAdam and E. van Sebille, “Surface connectivity and intero-cean exchanges from drifter-based transition matrices,” Journalof Geophysical Research: Oceans , 514–532 (2018). A. Gordon and R. Fine, “Pathways of water between the Pacificand Indian oceans in the Indonesian seas,” Nature , 46–149(1996). D. Tillinger and A. L. Gordon, “Fifty Years of the In-donesian Throughflow,” Journal of Climate , 6342–6355 (2009), https://journals.ametsoc.org/jcli/article-pdf/22/23/6342/3951874/2009jcli2981_1.pdf. R. Tarjan, “Depth-first search and linear graph algorithms,”SIAM J. Comput. , 146–160 (1972). M. Kubota, “A mechanism for the accumulation of floating ma-rine debris north of Hawaii,” J. Phys. Oceanogr. , 1059–1064(1994). 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. Bram-bini, and J. Reisser, “Evidence that the great pacific garbagepatch is rapidly accumulating plastic,” Scientific Reports , 4666(2018). J. R. Jambeck, R. Geyer, C. Wilcox, T. R. Siegler, M. Perry-man, A. Andrady, R. Narayan, and K. L. Law, “Plastic wasteinputs from land into the ocean,” Science , 768–771 (2015),https://science.sciencemag.org/content/347/6223/768.full.pdf. R. Yamashita and A. Tanimura, “Floating plastic in the KuroshioCurrent area, western North Pacific Ocean,” Marine PollutionBulletin , 485–488 (2007). M. van der Mheen, C. Pattiaratchi, and E. van Sebille, “Roleof indian ocean dynamics on accumulation of buoyant debris,”Journal of Geophysical Research: Oceans , 2571–2590 (2019). 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 particledynamics,” Geophys. Res. Lett. in press (2020). A. Cozar, E. Marti, C. M. Duarte, J. Garcia-de Lomas, E. vanSebille, T. J. Ballatore, V. M. Eguiluz, J. I. Gonzalez-Gordillo,M. L. Pedrotti, F. Echevarria, R. Trouble, and X. Irigoien, “Thearctic ocean as a dead end for floating plastics in the north at-lantic branch of the thermohaline circulation,” Science Advances , e1600582 (2017). A. Mountford and M. A. Morales Maqueda, “Eulerian modellingof the three-dimensional distribution of seven popular microplas-tic types in the global ocean,” Journal of Geophysical Research:Oceans , 1–16 (2019). 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 biogeo-graphic classification of the mesopelagic zone,” Deep Sea Re-search , 85 – 102 (2017). P. Ryan, “A simple technique for counting marine debris at seareveals steep litter gradients between the Straits of Malacca andthe Bay of Bengal,” Marine pollution bulletin , 128–136 (2013). M. Dellnitz, G. Froyland, C. Horenkam, K. Padberg-Gehle, andA. Sen Gupta, “Seasonal variability of the subpolar gyres in thesouthern ocean: a numerical investigation based on transfer op-erators,” Nonlinear Process. Geophys. , 655–663 (2009). FIG. 10. Reactive currents from each subtropical gyre patchto all other subtropical gyre patches. The source set is indi-cated in black in each panel; the target sets are indicated inred. Note the difference in the scales across the panels. FIG. S1. As in Fig. 4, but assuming α = in (28).FIG. S2. As in Fig. 4, but assuming α = in (28).FIG. S3. As in Fig. 4, but assuming α = 1 in (28). FIG. S4. As in Fig. 5, but assuming α = in (28).FIG. S5. As in Fig. 5, but assuming α = in (28).FIG. S6. As in Fig. 5, but assuming α = 1= 1