Dynamical geography and transition paths of Sargassum in the tropical Atlantic Open Access

Since early 2010s, pelagic Sargassum, a genus of brown macroalgae that form floating rafts at the ocean surface, has inundated the Intra-American Seas (IAS), particularly the Caribbean Sea, during the spring and summer months (Wang et al., 2019). These rafts of algae serve as important habitats for diverse marine fauna (Bertola et al., 2020) and have been argued to be important carbon sinks (Paraguay-Delgado et al., 2020). Nonetheless, these rafts carry large amounts of toxic substances and heavy metals, and, when they enter coastal zones, they can result in mortality to fishes and sea turtles and effectively smother seagrass and coral communities (van Tussenbroek et al., 2017). As large Sargassum rafts decompose onshore, they give off toxic amounts of hydrogen sulfide that can cause health problems in humans, the volume that washes ashore on beaches also diminishes tourism and, as a result, disrupts the local economy (Smetacek and Zingone, 2013; Resiere et al., 2018).

Comprehensive analyses of satellite imagery across the Atlantic Ocean have revealed a recurrent band of Sargassum high density between about 5 and 10°N, referred to as the Great Atlantic Sargassum Belt (GASB) in Wang et al. (2019), often extending off the coast of West Africa to the Gulf of Mexico (Fig. 1, bottom panel). The factors that precipitated blooms of pelagic Sargassum in the tropical Atlantic remain an area of active debate, as do the factors that maintain its occurrence across this region (Jouanno et al., 2021b; 2021a; and Lapointe et al., 2021). Several published studies have attempted to account for the spatiotemporal variability in its distribution by simulating the movement of Sargassum rafts as passive tracers advected by surface velocity fields produced by ocean circulation models. Conclusions about the connectivity between Sargassum blooms in the tropical Atlantic and the Sargasso Sea fundamentally differ depending on whether simulations assume winds contribute to Sargassum movement (Franks et al., 2016; Wang et al., 2019; and Johns et al., 2020). Available empirical data suggest that accounting for windage effects improves predictions of Sargassum raft trajectories (Putman et al., 2020) and distribution (Berline et al., 2020).

The importance of windage in modeling Sargassum movement is demonstrated in the top and middle panels of Fig. 1, which show discrete probability densities of finding Sargassum in the tropical Atlantic after three months of evolution from an initially uniform distribution inside $I=$ [5°S, 18°N] × [58°W, 15°W] as in Wang et al. (2019). Unlike Wang et al. (2019) and earlier work, which resorted to tracer advection by simulated velocity, the evolution is here provided by a transfer operator (Lasota and Mackey, 1994) obtained by discretizing (Ulam, 1960) the motion as described by the trajectories x(t) of surface drifting buoys from the National Oceanic and Atmospheric Administration (NOAA) Global Drifter Program (GDP) (Lumpkin and Pazos, 2007) assuming that the underlying dynamics represent a passive advection–diffusion process [e.g. Miron et al. (2019a)].

The discrete transfer operator is given by a transition matrix $P=(Pij)∈RN×N$ of conditional probabilities of x(t) to moving between nonintersecting boxes {b₁, …, b_N} covering the surface of the Atlantic between 9°S and 39°N, denoted $D$⁠, viz.,

where

independent of the start time t. Here, X_n denotes random position at discrete time nT, where T is the transition time step. The boxes represent the states of a Markov chain, namely, a stochastic model describing the stochastic state transitions in which the transition probability of each state depends only on the state attained in the previous event (Norris, 1998). By ignoring the start time t of a trajectory in (1), the chain is rendered time-homogeneous. The evolution of prob(X_n ∈ b_i) is obtained under left multiplication by $P$⁠, i.e.,

Figure 1 specifically shows prob(X₂₄ ∈ b_i), where prob(X₀ ∈ b_i) = 1/|I| for b_i∈I spanning $I⊂D$ and 0 otherwise. Assuming T = 5 d, as we do here, this represents the end of nearly four months of evolution. Our transition time step choice T = 5 d is long enough for Markovian dynamics to approximately hold, as it is longer than the typical Lagrangian integral time scale, which is about 1 day at the ocean surface, e.g. LaCasce (2008), and as noted in pioneering (Maximenko et al., 2012) and most recent work [Miron et al. (2017), Olascoaga et al. (2018), McAdam and van Sebille (2018), Miron et al. (2019b;,2019a; 2021), and Drouin et al. (2022)], involving observed, rather than simulated, trajectory data.

A critical difference between the top and middle panels of Fig. 1 is that the former uses trajectories of drifters drogued at 15 m, while the latter uses trajectories of undrogued drifters (Lumpkin et al., 2012). Note that the assessment based on undrogued drifter motion (middle panel) much better resembles the GASB (bottom panel) than the assessment based on drogued drifter motion (top panel). Similar behavior has been recently reported by van Sebille et al. (2021) from direct inspection of individual trajectories produced by a reduced set of custom-made undrogued drifters along with drogued drifters of the GDP type. As opposed to drogued GDP drifter motion, which mostly reflects water motion at 15 m depth in the water column, undrogued drifter motion is affected by ocean current shear between 0 and 15 m, and wind and wave drag effects, mediated by “inertial” effects, i.e., those due to the buoyancy, size, and shape of the drifters (Furnans et al., 2008; Beron-Vera et al., 2016; 2019; Olascoaga et al., 2020; Miron et al., 2020b; 2020a; and Beron-Vera, 2021). The resemblance—impressively quite close—between the observed GASB distribution (in the bottom panel of Fig. 1) and that suggested by undrogued drifter trajectories indicates a dominant role of inertial effects in the drift of Sargassum. Indeed, these appear as an important mechanism for pelagic Sargassum to inundate the IAS.

The rest of this note is devoted to go beyond the discussion of direct probability density evolution with the twofold goal in mind to gain insight into: Sargassum connectivity as inferred by a Markov chain constructed using undrogued drifter trajectories exclusively, which show the best agreement with actual Sargassum movement (Sec. II); and how Sargassum pathways connecting a potential remote source with the IAS are accomplished in the most effective manner (Sec. III). This is done by applying specialized probabilistic tools, which are briefly reviewed before the results are discussed in each case. Section IV comments on aspects of the movement of Sargassum rafts that are needed to be accounted for in order to make the most out of the learnings from the analysis presented here. This brief communication is closed with a summary (Sec. V).

Time-asymptotic aspects of the evolution along the Markov chain resulting from the discretization of the drifter motion are encoded in the spectral properties of the transition matrix $P$⁠. If every state of the chain (box of the covering) is visited irrespective of the starting state [property called irreducibility or ergodicity, and formally represented as $∀i,j:bi,bj∈D,∃nij∈Z0+\{∞}:(Pnij)ij>0$] and no state is revisited cyclically [property referred to as aperiodicity or mixing and expressed as $∃i:bi∈D:gcd{n∈Z0+:(Pn)ii>0}=1$], then the eigenvalue λ = 1 of $P$ is the largest possible and has multiplicity one.

The left eigenvector, π, corresponding to λ = 1, is strictly (componentwise) positive and represents an invariant distribution (i.e., $π=πP$⁠), which typically reveals local maxima (bumps) where trajectories settle on in the long run. Normalizing π to a probability vector (i.e., such that $∑i:bi∈Dπi=1$⁠), it follows that $π=pP∞$ for any probability vector p. In other words, π is a limiting or stationary distribution that time-asymptotically sets prob(X_n ∈ b_i) = π_i.

The right eigenvector corresponding to λ = 1 is a (column) vector of ones, denoted 1. Indeed $P1=1$ due to (row) stochasticity of $P$⁠, viz., $∑j:bj∈DPij=1$ for all i, trivially indicating a backward-time basin of attraction for the time-asymptotic attractors where π maximizes: any probability vector with support in $D$ distributes in the long run like π.

Left eigenvectors of P with λ close to 1 are signed vectors that decay at the slowest rate. The regions where $P$ these locally maximize thus indicate almost-invariant attracting sets. Plateaus in the corresponding right eigenvectors highlight basins of attraction or regions where trajectories temporarily converging in the almost-invariant attractors begin. This imposes constraints on connectivity (Dellnitz and Hohmann, 1997; Dellnitz and Junge, 1999; Froyland and Dellnitz, 2003; and Koltai, 2010). The collection of nonoverlapping basins of attraction forms a partition of the dynamics into the weakly (dis)connected “provinces” of what has been termed (Froyland et al., 2014b; Miron et al., 2017) a dynamical geography, i.e., one that does not depend on arbitrary geographic demarcations but rather on the intrinsic flow characteristics of each province.

The eigenvector method, above, has been applied to surface drifter and submerged float trajectories suggesting a characterization of preferred pollution centers in the subtropical gyres (Froyland et al., 2014b; Miron et al., 2019b; and 2021) and in marginal seas, both at the surface (Miron et al., 2017) and at depth (Miron et al., 2019a), as almost-invariant attracting sets with corresponding basins of attraction spanning areas as large as those of the corresponding geographic basins, suggesting a strong influence of the regions collecting pollutants on their transport.

A similar strong influence imposed by the northwestern Gulf of Mexico, the Sargasso Sea, and the Gulf of Guinea is revealed from the analysis of our $P$⁠, constructed using undrogued GDP drifter trajectories from the NOAA dataset. The top panel of Fig. 2 shows the locations where the first and second subdominant left eigenvectors of P maximize, which identify the aforementioned geographical locations as almost-invariant attracting sets. The corresponding basins of attraction, as revealed by the first and second subdominant right eigenvectors of $P$⁠, are shown in the bottom panel of Fig. 2. These are depicted with the same color as the corresponding attractors to facilitate the association. [In every case, we have made use of the sparse eigenbasis approximation (Froyland et al., 2019), which enables multiple feature extraction with application to coherent set identification.] Note the large geographical regions covered by the basins, which represent the domains of influence of the corresponding attractors. For instance, trajectories passing through the Caribbean Sea on their way into the northwestern Gulf of Mexico most likely start anywhere in the Atlantic domain shown here except in the subtropical gyre. This suggests weak communication between the Sargasso Sea and the GASB region, at least for the average undrogued GDP drifter motion.

To unveil specific connectivity routes, additional probabilistic tools are needed. Particularly useful are those provided by the Transition Path Theory (TPT) (E and Vanden-Eijnden, 2006; Metzner et al., 2009), which has been recently adapted to open dynamical systems (Miron et al., 2021).

Developed to investigate rare events in complex systems, such as chemical reactions or conformation changes of molecules, TPT provides a statistical characterization of ensembles of “reactive” trajectories, namely, trajectories along which direct transitions between a source set (or “reactant”) A and a target set (or “product”) B in a Markov chain take place (Fig. 3).

Applications of TPT have now gone beyond the study of molecular dynamics or chemical kinetics (Noé et al., 2009; Metzner et al., 2006; 2009; Meng et al., 2016; Thiede et al., 2019; Liu et al., 2019; and Strahan et al., 2021). TPT has been used to shed light on pollution pathways in the ocean (Miron et al., 2021), paths of the upper (Drouin et al., 2022) and lower (Miron et al., 2022) limbs of the overturning circulation in the Atlantic Ocean, and even atmospheric phenomena such as sudden stratospheric warmings (Finkel et al., 2020).

The main objects of TPT are the forward, $q+=(qi+)∈R1×N$⁠, and backward, $q−=(qi−)∈R1×N$⁠, committor probabilities. These give the probability of a random walker initially in b_i to first enter B and last exit A, respectively. The committors are fully computable from $P$ and π, according to

Here, $P+=P$⁠; $Pij−≔prob(Xn∈bj∣Xn+1∈bi)=πjπiPji$ are the entries of the time-reversed transition matrix, i.e., for the original chain {X_n} traversed in backward time, {X_−n}; and |_S denotes restriction to the subset S while |_S,S′ that to the rows corresponding to S and columns to S′.

The committor probabilities are used to express several statistics of the ensemble of reactive trajectories. The main ones are as follows:

1. The distribution of reactive trajectories, $πAB=(πiAB)∈R1×N$⁠, where $πiAB$ is defined as the joint probability that a trajectory is in box b_i while transitioning from A to B, and is computable as [Metzner et al. (2009) and Helfmann et al. (2020)]

   $πiAB=qi−πiqi+.$

   (5)

   It describes the bottlenecks during the transitions, i.e., where reactive trajectories spend most of their time. Clearly, $πi:bi∈A∪BAB≡0$⁠.
2. The effective current of reactive trajectories, $f+=(fij+)∈RN×N$⁠, where $fij+$ gives the net flux of trajectories going through b_i at time nT and b_j at time (n + 1)T on their way from A to B, indicates the most likely transition channels. This is computable according to (Helfmann et al., 2020)

   $fij+=maxfijAB−fjiAB,0,fijAB=qi−πiPijqj+.$

   (6)

In the top panel of Fig. 4, we show the effective current of reactive trajectories resulting between a hypothetical source set A off West Africa (blue box) and a target set B formed by the union of boxes covering the Gulf of Mexico (red boxes) when TPT is applied on the Markov chain constructed using undrogued drifter trajectories. To visualize $f+$⁠, we follow Helfmann et al. (2020) and Miron et al. (2021) and for each b_i we estimate the vector of the average direction and the amount of the effective reactive current to each b_j, j ≠ i. The location of A has been intentionally chosen to lie in a region near the triple boundary of the dynamical provinces (cf. Fig. 2, bottom panel), which minimizes the influence of the almost-invariant forward attractors within each basin of attraction: a probability density initialized there, in principle, has equal chance to converge, temporarily, into any of each of the three almost-attracting sets. On the other hand, placing A off West Africa as chosen finds rationale in the fact that upwelling-favorable winds there may provide the required nutrients to trigger blooming by vertical pumping them into the mixed layer. Note that the eastern extent of the GASB coincides roughly with this location (cf. Fig. 1). The bottom panel of Fig. 4 uncovers the spatial locations of the bottleneck of the transition paths, namely, where the paths spend most of their time in their direct transition from the source into the target.

Note the two routes revealed in the top panel of Fig. 4: a direct route roughly along the GASB and another more southern route, reported for the first time here, going through the Gulf of Guinea, and then passing by the mouth of the Amazon River. The southern path is strongly influenced by the almost-invariant attracting set in the Gulf of Guinea, which makes the trajectories bottleneck there. The latter is highlighted by the large values taken by the distribution of reactive trajectories in the bottom panel of Fig. 4.

Two additional TPT diagnostics are as follows:

1. The rate of reactive trajectories entering B (or, equivalently, exiting A), defined as the probability per time step of reactive trajectory to enter B (or exit A), which is computed as (Metzner et al., 2009 and Helfmann et al., 2020)

   $kAB=∑i,j:bi,bj∈BfijAB.$

   (7)

   Dividing by the transition time step T, k^AB has the interpretation of the frequency at which a trajectory enters B (or exists A) (Miron et al., 2021).
2. The expected duration of a transition from A to B, which is obtained by dividing the probability of a trajectory piece being reactive by the transition rate (interpreted as a frequency) (Helfmann et al., 2020),

   $tAB=∑j:bj∈D\A∪BμjABkAB.$

   (8)

The computed expected duration of the transition paths into the IAS is about 12 yr, which seems rather long. However, this should not come as big surprise, given the strong influence exerted by the Gulf of Guinea almost-attracting set. In fact, if the transitions are set to avoid the basin of attraction of Gulf of Guinea attractor, the expected duration drops down to 1.8 yr or so.

The reactive currents in the case that the basin of attraction of the Gulf of Guinea attractor is avoided reveal a single route into the IAS, roughly aligned along the GASB. The physical significance of this almost-invariant attractor is demonstrated in the top panel Fig. 5, which provides proof, independent of the Markov chain model, of the existence of such an attractor. The figure shows individual trajectories of 64 undrogued GDP drifters from the NOAA database over the period 1998–2014, which were found to converge in the Gulf of Guinea. A previous TPT analysis (Miron et al., 2021) had already identified the neighborhood of the Gulf of Guinea as a province for a committor-based dynamical geography for global marine debris circulation. It was also pointed out (Sutton et al., 2017) that the Gulf of Guinea is a mesopelagic niche with genomic characteristics that are different than its surroundings.

The southern GASB branch unveiled by the TPT analysis does not appear to be an artifact of it. Support for the southern transition path is given in the bottom panel Fig. 5, which shows the fourth-root transformation of the quantity shown in the bottom panel of Fig. 1 computed from boreal summer composites over 2011–2018. Note the Sargassum presence below the equatorial line, albeit in less abundance than along the (main) GASB (branch). This is consistent with the TPT results (Fig. 4), which reveals that the southern path is weaker than the northern path, where the probability of the reactive currents is stronger consistent with the higher concentration of satellite-inferred floating Sargassum.

The analysis of the Markov chain model derived using undrogued drifter trajectories indicates that the combined effects of winds, waves, and surface currents can account for the basin-scale distribution of the GASB. These results are in agreement with several other studies that have suggested the importance of windage and currents in shaping the trajectories and distribution of Sargassum (Brooks et al., 2018; Putman et al., 2020; and Berline et al., 2020). Our probabilistic analysis also indicates three, mostly independent, oceanographic provinces associated with aggregations of Sargassum that occur in the Sargasso Sea, Gulf of Guinea, and Gulf of Mexico; a finding consistent with Wang et al. (2019). Furthermore, our analysis suggests the existence of a secondary pathway by which Sargassum is transported from the tropical Atlantic into the Caribbean Sea, south of the GASB (Figs. 4 and 5). While this pathway has not been previously described, evidence consistent with this pathway is suggested in earlier studies. Notably, Putman et al. (2018) showed the importance of the Guinea Current for moving Sargassum along the coast of northern Brazil into the Caribbean, which is fed by the South Equatorial Current. Satellite imagery depicted in Putman et al. (2018) and Gower et al. (2013) suggests that Sargassum enters this region from south of the equator. Thus, investigation of the ocean–atmosphere dynamics associated with this newly described southern, slower pathway for Sargassum and the more northern, faster GASB (Wang et al., 2019), could improve the ability to predict basin-scale changes in Sargassum distribution, especially given very different time scales associated with the two pathways.

For instance, it may be important to consider that the North Equatorial Counter Current, running between 5 and 10°N, exhibits seasonal reversal (Lumpkin and Garzoli, 2005), which is not taken into account by our autonomous TPT analysis, wherein the transition matrix construction ignores the starting day of a drifter trajectory. To investigate the influence of the seasonal variability in our results, a nonautonomous TPT analysis will be needed. The theoretical basis for such an analysis is now available (Helfmann et al., 2020), but still to be implemented. A main constraint in this case of interest is imposed by the data availability, which is unlikely to resist a seasonal partition. Simulated trajectories may serve the purpose, but this is beyond our scope here as is assessing interannual variability in Sargassum connectivity,

Beyond the assumed time-homogeneity, the results of this work are limited in two important counts. First is the physiological changes that Sargassum rafts experience as they are dragged by the combined action of ocean currents and winds, which the Markov chain model constructed here did not account for. A second aspect is the architecture of the Sargassum rafts, which plays a role in advection. It is well-known that finite-size or “inertial” particles, even neutrally buoyant ones, immersed in a fluid cannot adapt their velocities to the carrying flow (Maxey and Riley, 1983). Floating particle motion is further complicated by the fact that the carrying flow is not given by the ocean velocity below the sea surface itself, but by the latter plus a buoyancy-dependent fraction of the wind velocity right above the sea surface (Beron-Vera et al., 2019; Olascoaga et al., 2020; Miron et al., 2020b; 2020a; and Beron-Vera, 2021). This aspect is implicitly accounted for by our Markov chain model as the undrogued drifters are affected by inertial effects. However, Sargassum rafts are not isolated inertial particles. Rather, they can be more precisely envisioned as elastic networks of buoyant inertial particles, with the buoyant inertial particles representing the gas-filled bladders connected by flexible stems that keep the rafts afloat. This was not represented by our Markov chain model. A mathematical model has been recently proposed by Beron-Vera and Miron (2020) and used Andrade-Canto et al., 2022) to explain observations. This model, coupled with a physiological model of the transformations that a Sargassum raft experiences as it travels across regions with varying environmental conditions, would likely further elucidate the complex dynamics of Sargassum invasions in the IAS. Such work could further improve transport predictions by distinguishing between the roles of physical and physiological processes on observed distributions.

In this note, we have constructed a time-homogeneous Markov chain from trajectories of undrogued drifters from the NOAA Global Drifter Program. With this tool, we found that probability density can distribute in a manner that resembled very closely the recently observed recurrent belt of high Sargassum concentration in the tropical Atlantic that has been termed the Great Atlantic Sargassum Belt (GASB). A spectral analysis of the associated transition matrix was further found to unveil a forward attracting almost-invariant set in the northwestern Gulf of Mexico with a corresponding basin of attraction weakly connected with the Sargasso Sea, but including the nutrient-rich regions at the outflow of the Amazon and Orinoco Rivers and also, marginally, the upwelling system off the northern coast of west Africa. This represents a data-based inference of potential remote sources of Sargassum recurrently invading the Intra-Americas Seas (IAS). By further applying Transition Path Theory (TPT) on the data-derived Markov chain model, two potential pathways for Sargassum into the IAS from the upwelling system off the coast of Africa were revealed. One TPT-inferred pathway was found to take place along the GASB. The second pathway, much slower than the GASB, was found to be more southern, first going through the Gulf of Guinea, then across the Atlantic toward the mouth of the Amazon River, and finally along the eastern South American margin. Why satellite imagery does not reveal an intense southern GASB branch needs to be investigated in view of its potentially important consequences for bloom stimulation by nutrients from river runoff. This should be done by accounting for the physiological transformations that Sargassum rafts are acted upon as they are transported by the combined action of currents and winds, which is the subject of ongoing research.

We are grateful to P. Miron for the benefit of many helpful discussions relating time-asymptotic almost-invariant sets and transition path theory. This work was supported by the University of Miami’s Cooperative Institute for Marine and Atmospheric Studies and National Science Foundation Grant No. OCE2148499.

The authors declare no competing interest.

F. J. Beron-Vera: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). M. J. Olascoaga: Conceptualization (supporting); Formal analysis (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). N. F. Putman: Conceptualization (equal); Writing – review & editing (equal). J. Triñanes: Data curation (supporting). G. J. Goni: Funding acquisition (equal). R. Lumpkin: Funding acquisition (equal).

The data that support the findings of this study are openly available in NOAA Global Drifter Program at https://www.aoml.noaa.gov/phod/gdp/index.php (Lumpkin and Pazos, 2007) and USF and distributed by SaWS at https://optics.marine.usf.edu/projects/saws.html (Wang and Hu, 2016).