Charting the course of Sargassum: Incorporating nonlinear elastic interactions and life cycles in the Maxey–Riley model

Abstract

The surge of pelagic Sargassum in the Intra-America Seas, particularly the Caribbean Sea, since the early 2010s has raised significant ecological concerns. This study emphasizes the need for a mechanistic understanding of Sargassum dynamics to elucidate the ecological impacts and uncertainties associated with blooms. By introducing a novel transport model, physical components such as ocean currents and winds are integrated with biological aspects affecting the Sargassum life cycle, including reproduction, grounded in an enhanced Maxey–Riley theory for floating particles. Nonlinear elastic forces among the particles are included to simulate interactions within and among Sargassum rafts. This promotes aggregation, consistent with observations, within oceanic eddies, which facilitate their transport. This cannot be achieved by the so-called leeway approach to transport, which forms the basis of current Sargassum modeling. Using satellite-derived data, the model is validated, outperforming the leeway model. Publicly accessible codes are provided to support further research and ecosystem management efforts. This comprehensive approach is expected to improve predictive capabilities and management strategies regarding Sargassum dynamics in affected regions, thus contributing to a deeper understanding of marine ecosystem dynamics and resilience.

Significance Statement

The massive increase in Sargassum seaweed in the Caribbean Sea since the early 2010s poses major ecological challenges. This research introduces a new model that integrates ocean currents, winds, and biological factors to better understand and predict Sargassum movement. Unlike existing models, this new approach simulates interactions within Sargassum rafts and their aggregation in ocean eddies, providing a more accurate prediction of their movement. Validated with satellite data, the model surpasses current methods, offering improved tools for researchers and policymakers. The model is distributed in a fully open source manner and features an interface that enables use by nonexperts.

“Just before it was dark, as they passed a great island of Sargasso weed that heaved and swung in the light sea as though the ocean were making love with something under a yellow blanket, ….”

Ernest Hemingway, The Old Man and the Sea.

Introduction

The brown macroalgae Sargassum was reportedly first documented by Christopher Columbus and his crew in 1492 (1). The floating rafts, characterized by small gas-filled bladders, apparently reminded them of “salgazo,” a grape variety found in Portugal. Originally, abundant throughout the North Atlantic subtropical gyre, known as the Sargasso Sea. Subsequently, beginning in the early 2010s, massive pelagic Sargassum rafts began to inundate the Intra-America Seas, notably the Caribbean Sea, during the spring and summer months (2).

These pelagic Sargassum rafts provide essential habitats for a diverse range of marine fauna, significantly enhancing the region’s biodiversity (3). Additionally, they act as carbon sinks, contributing to the global carbon cycle (4). However, these benefits come with notable drawbacks. The rafts accumulate toxic substances and heavy metals. Upon entering coastal zones, they can cause the mortality of fish and sea turtles and suffocate vital seagrass and coral communities (5). The breakdown of extensive onshore Sargassum rafts emits hydrogen sulfide, which may present health hazards to people. Furthermore, the mass stranding of algae on beaches significantly reduces tourism, disrupting the local economy (6, 7).

Comprehensive analyses of satellite imagery have revealed a persistent band of high-density Sargassum across the tropical North Atlantic (8, 9). The factors that precipitate blooms of pelagic Sargassum remain a subject of active scientific debate and inquiry, as do the factors that maintain its presence across this region (9–11). Currently, the most significant source of uncertainty is the fundamental understanding of the movement of pelagic Sargassum (12–15). Without a mechanistic understanding of Sargassum dynamics, it is challenging to assess how much of the observed distributions result from transport processes versus local growth, and how these distributions will evolve over time (15–17).

Several studies have attempted to account for the spatiotemporal variability in its distribution by tracking synthetic fluid particles within ocean circulation model velocity fields. Results vary significantly depending on whether simulations assume winds contribute to Sargassum movement (e.g. through a windage or Stokes drift factor) and on assumptions about growth and mortality (typically accounted for in models by how long particles are tracked). Some studies indicate evidence of occasional transport from the Sargasso Sea into the tropics, while others suggest this to be unlikely (9, 10, 18).

To address the uncertainties associated with the largely piecemeal and ad hoc approaches in Sargassum modeling so far, this article develops a mechanistic model that integrates two critical aspects: one physical and one biological. The physical aspect involves transport resulting from the combined action of ocean currents and winds, mediated by inertial effects. The biological aspect encompasses Sargassum reproduction through vegetative growth and fragmentation and death through factors such as aging and sinking, mediated by general mortality.

The modeling of inertial effects builds on the Maxey–Riley theory for the motion of finite-size, buoyant (i.e. inertial) particles immersed in fluid flow (19–21). The Maxey–Riley equation represents an application of Newton’s second law, incorporating various forces, including the flow force exerted on the particle by the undisturbed fluid, the added mass force resulting from part of the fluid moving with the particle, and the drag force caused by the fluid viscosity, among others. These forces are derived from first principles, as they originate by integrating the Cauchy stress tensor over the spherical surface of the particle immersed in a moving viscous fluid, assuming the Stokes number is small (21). For a turbulent flow such as the ocean flow, this imposes a restriction on the particle size; depending on the chosen scales, its radius should be (much) smaller than about 1 meter.

The motion of inertial particles thus deviates fundamentally from that of infinitesimally small, neutrally buoyant (i.e. fluid) particles, as inertial particles are incapable of adjusting their velocities to conform to the flow velocity (22). Although the trajectories of neutrally buoyant particles eventually synchronize with those of fluid particles over time, they do so relative to fluid particles starting from different initial positions (23, 24). The manifestation of inertial effects, which have been documented even for minute neutrally buoyant particles (25), can be expected to be significantly magnified in the dynamic context of sprawling rafts of floating Sargassum.

The fluid mechanics Maxey–Riley theory has been extended to particles floating at the ocean surface (26), enabling them to interact through elastic forces (27). The oceanographic Maxey–Riley framework for individual particles has undergone successful testing in both the field (28, 29) and the laboratory (30) and is exhaustively reviewed in (31).

Elastically interacting particles serve as a parsimonious model for Sargassum plant architecture, representing flexible stems as springs and air-filled bladders as inertial particles. These interactions lead to behavior that differs from noninteracting floating inertial particles, promoting concentration within mesoscale eddies independent of their polarity. This behavior is supported by rigorous results (27) and empirical evidence (32).

The Maxey–Riley model for Sargassum plant motion is significantly extended here to simulate interactions among rafts, conceptualized as networks subjected to nonlinear elastic forces. The modeling of Sargassum reproduction is approached mechanistically by regulating the size of the networks based on a specified law.

In a notable stride forward, this article not only validates the model against satellite-derived Sargassum densities, outperforming the so-called leeway approach to transport that supports contemporary Sargassum modeling (14, 33, 34). Additionally, this article develops and publicly releases codes for executing model simulations and computing densities. Consequently, the community will have access to a robust and flexible toolkit for in-depth exploration of Sargassum dynamics, enhancing understanding and predictive capabilities in Sargassum research, management, and response efforts.

Results

Model overview

We delineate the terminology pertinent to directly related models. The BOM model was introduced in (26) as a comprehensive Maxey–Riley framework to characterize the dynamics of floating particles of small but finite dimensions on the ocean surface. The eBOM model was presented in (27) as a model for Sargassum plants, conceptualized as networks of elastically interacting BOM particles. In this article, we advance the eBOMB model as an extension of the aforementioned model, as elaborated in the subsequent sections.

The fundamental object of the eBOMB model is a clump of Sargassum in the open ocean, subject to currents and wind. A clump is qualitatively a handful of Sargassum that we consider to be a hard sphere with a small but finite radius. We call a network of interacting clumps a raft and assume that the motion of a continuous macroscopic quantity of Sargassum can be well approximated by the motion of this discrete set of clumps as “landmarks.” This is represented schematically along with a comparison to the precursor eBOM model in Fig. 1.

Fig. 1.

A, B) As viewed from above, an individual Sargassum plant is idealized as a network of air bladders connected by semi-rigid stipes. The eBOM model treats this as a network of small inertial spheres, each of which obeys a Maxey–Riley equation coupled by linear elastic (Hookian) forces. C, D) As viewed from above, a Sargassum raft is idealized as a continuous biomass dotted with indicator clumps that approximate that raft’s motion. The eBOMB model treats the clumps as Maxey–Riley particles similarly to eBOM, but the connecting springs have nonlinear stiffness parameters that allow for fragmentation and recombination of subrafts.

The totality of Sargassum under consideration at a given time is itself a raft although it may consist of a number of disconnected networks of clumps which may be referred to as subrafts. The overall size hierarchy of Sargassum quantities is therefore, from smallest to largest: air bladders, stipes, plants, clumps, subrafts, and rafts. The model consists of three main components, the clump physics controlled by coupled Maxey–Riley equations, spring forces with nonlinear stiffness parameters connecting subrafts, and the biological model controlling the growth and death of clumps.

Physical module of the eBOMB model

The Maxey–Riley physics follows (26, 27). There are four parameters required to define the main equations: α is the unitless parameter measuring windage, τ is the inertial response time, and R is a dimensionless parameter that quantifies the exposure of the assumed spherical clump to the air. These parameters depend further on the radius a of a clump and the density ratio of the clump to the water $\delta \ge 1$⁠, which we refer to as buoyancy (additional details are given in Materials and methods). We will assume here that each of these parameters can be specified independently, although we will justify their particular values with reference to δ and a later. Let $\mathbf{x}=(x,y)$ denote the position on the ocean surface with rescaled longitude–latitude (i.e. geographic) coordinates $x=a_{\odot }\text{cos}{\vartheta }_0(\lambda -{\lambda }_0)$ and $y=a_{\odot }(\vartheta -{\vartheta }_0)$⁠, where $a_{\odot }$ is the mean Earth’s radius and $({\lambda }_0,{\vartheta }_0)$ denotes a reference location. We define ${\gamma }_{\odot }(y)=\text{sec}{\vartheta }_0\text{cos}\vartheta$ and ${\tau }_{\odot }(y)=a_{\odot }^{-1}\tan \vartheta$⁠, which are geometric factors that arise due to the sphericity of the Earth. Let $\mathbf{v}(\mathbf{x},t)$ and ${\mathbf{v}}_{\mathrm{W}}(\mathbf{x},t)$ denote the near-surface ocean and wind velocities, respectively. Define

$$\mathbf{u}=(1-\alpha )\mathbf{v}+\alpha {\mathbf{v}}_{\mathrm{W}}.$$

(1)

We decompose $\mathbf{v}={\mathbf{v}}_{\mathrm{E}}+\sigma {\mathbf{v}}_{\mathrm{S}}$⁠, where ${\mathbf{v}}_{\mathrm{E}}$ represents the part of the ocean velocity, sometimes referred as “Eulerian,” on time scales much longer than the period of surface gravity waves and ${\mathbf{v}}_{\mathrm{S}}$ is the wave-induced Stokes drift. This decomposition follows (35), except for the unitless quantity σ, which we introduce to parameterize uncertainty around ${\mathbf{v}}_{\mathrm{S}}$⁠. Let $\omega ={\gamma }_{\odot }^{-1}{\mathrm{\partial }}_x{v}_y-{\mathrm{\partial }}_y{v}_x+{\tau }_{\odot }{v}_x$ be the vorticity of the ocean velocity and let ${\mathrm{\nabla }}_{\mathbf{w}}\mathbf{w}={\gamma }_{\odot }^{-1}{w}_x{\mathrm{\partial }}_x\mathbf{w}+w_y{\mathrm{\partial }}_y\mathbf{w}+{\tau }_{\odot }{w}_x{\mathbf{w}}^{\perp }$ be the covariant derivative of $\mathbf{w}=\mathbf{v}$ or $\mathbf{u}$ in the direction of itself, where ⊥ stands for anticlockwise $\frac{\pi }{2}$-rotation. The trajectory ${\mathbf{x}}_i(t)$ of the ith clump of a raft obeys (see Supplementary material for details):

$$\sqrt{\mathsf{m}}{\dot{\mathbf{x}}}_i=\mathbf{u}{|}_i+\tau {\mathbf{u}}_{\tau }{|}_i+\tau {\mathbf{F}}_i,$$

(2a)

where

$${\mathbf{u}}_{\tau }=R({\mathrm{\partial }}_t\mathbf{v}+{\mathrm{\nabla }}_{\mathbf{v}}\mathbf{v})+R(f+\frac{1}{3}\omega ){\mathbf{v}}^{\perp }-{\mathrm{\partial }}_t\mathbf{u}-{\mathrm{\nabla }}_{\mathbf{u}}\mathbf{u}-(f+{\tau }_{\odot }{u}_x+\frac{1}{3}R\omega ){\mathbf{u}}^{\perp }.$$

(2b)

In this context, the matrix $\mathsf{m}$ is the representation of metric in the (rescaled geographic) coordinates $(x,y)$⁠, with elements ${\mathsf{m}}_{11}={\gamma }_{\odot }^2$⁠, ${\mathsf{m}}_{11}=1$ and ${\mathsf{m}}_{12}=0={\mathsf{m}}_{21}$⁠. The notation ${|}_i$ signifies the association with clump i. The function $f(y)=2\mathrm{\Omega }\sin \vartheta$⁠, where Ω denotes the Earth’s rotation rate, is the Coriolis parameter. Additionally, ${\mathbf{F}}_i$ denotes an external force. Consequently, for a raft consisting of N clumps, we obtain a system of $2N$ first-order ordinary differential equations, which are coupled via  ${\mathbf{F}}_i$⁠.

We now describe the interaction spring forces. Let ${\mathbf{F}}_i$ be the force felt by the clump labeled i and let ${\mathbf{x}}_{ij}={\mathbf{x}}_i-{\mathbf{x}}_j$⁠, that is,

$${\mathbf{F}}_i=-{\displaystyle \sum _{j\in \mathrm{neighbor}(i,t)}}k(|{\mathbf{x}}_{ij}|)(|{\mathbf{x}}_{ij}|-L)\frac{{\mathbf{x}}_{ij}}{|{\mathbf{x}}_{ij}|},$$

(3)

where $|{\mathbf{x}}_{ij}{|}^2={\gamma }_{\odot }^2(x_i-x_j{)}^2+(y_i-y_j{)}^2$⁠. Additionally, $k(|{\mathbf{x}}_{ij}|)$ represents the stiffness function, L is the natural length of the spring, and $\text{neighbor}(i,t)$ is equal to the set of indices of clumps that are connected to i at time t. Note that all springs have the same k and L. The motivation behind this force is the effect of inter-clump entanglement which we assume provides a restorative force up to a certain distance, at which the clumps dissociate completely. Hence, we write

$$k(|{\mathbf{x}}_{ij}|)=\frac{A}{{\mathrm{e}}^{(|{\mathbf{x}}_{ij}|-2L)/\mathrm{\Delta }}+1},$$

(4)

where Δ is taken small enough such that $k(|{\mathbf{x}}_{ij}|)\approx A$ for $0\le |{\mathbf{x}}_{ij}|\le 2L$ and $k(|{\mathbf{x}}_{ij}|)\approx 0$ for $|{\mathbf{x}}_{ij}|>2L$⁠. It is worth mentioning that a piecewise constant $k(|{\mathbf{x}}_{ij}|)$ produces almost the same numerical results, yet Eq. 4 is a differentiable function that meets the technical conditions of Theorem 1. We take $\text{neighbor}(i,t)$ to be a set of nearest-neighbors to clump i at time t. This allows us to ensure that the simulations are fast by ignoring distant clumps, while also providing control over the size of subrafts, i.e. taking the number of nearest-neighbors to be large results in subrafts with greater membership. The natural length of the spring L is computed based on the initial configuration of clumps so that the model scales appropriately. In particular, if $d_i(K)$ is the mean distance of clump i’s K nearest-neighbors at the initial integration time, then L is chosen to be the median among all $d_i(K)$⁠.

Biological module of the eBOMB model

To simulate Sargassum life cycles, we employ a hybrid version of the models proposed by (12, 36), specifically adapted to the scenario of a discrete set of clumps as delineated below. Each clump is allocated an amount parameter, denoted in arbitrary units as $S(t)$⁠, and is initialized with the condition $S(0)=S_0$⁠. At every integration time step h of the coupled system defined by Eq. 2, $S(t)$ is updated via

$$S(t+h)=S(t)+(g(t)-m)h,$$

(5)

where $g(t)$ and m represent the factors controlling growth and mortality, respectively. We take m to be a constant and $g(t)$ equal to a constant multiplied by factors depending on the temperature T and nitrogen content of the water N via

$$g(t)={\mu }_{max}\cdot \frac{\mathcal{T}(t)}{k_{\mathrm{N}}/N(t)+1},$$

(6a)

where

$$\mathcal{T}(t)=\{\begin{array}{cc}\text{exp}(-{\displaystyle \frac{1}{2}{\left({\displaystyle \frac{T-T_0}{T-T_{min}}}\right)}^2})& \text{if}{T}_{min}\le T(t)\le T_0,\\ \text{exp}(-{\displaystyle \frac{1}{2}{\left({\displaystyle \frac{T-T_0}{T-T_{max}}}\right)}^2})& \text{if}{T}_0<T(t)\le T_{max},\\ 0\hfill & \text{otherwise}.\hfill \end{array}$$

(6b)

Here, ${\mu }_{max}$ is the maximum Sargassum growth rate, $T_{min},T_{max}$ are the minimum and maximum temperatures for Sargassum growth, $T_0=(T_{min}+T_{max})/2$ is the (approximate) optimal growth temperature and $k_{\mathrm{N}}$ is the Sargassum nitrogen uptake half saturation. Once Eq. 5 has been applied at the end of the time step, we compare $S(t+h)$ to thresholds $S_{min}$ and $S_{max}$⁠. If $S(t+h)>S_{max}$⁠, a new clump is born with $S_{\mathrm{new}}(t+h)=S_0$⁠. The new clump is placed a distance L away from its “parent” at an angle chosen uniformly at random, and the parent is reset to $S(t+h)=S_0$⁠. If $S(t+h)<S_{min}$⁠, the clump is removed completely from the integration by dynamically reducing the size of the integrator cache. Clumps may also die when reaching land.

Model validation

Our principal comparative framework employs the leeway model, which governs the transport dynamics of Sargassum within state-of-the-art simulations (14). Frequently referenced in the search-and-rescue literature (37), the leeway model is defined as

$${\dot{\mathbf{x}}}_i={\mathbf{v}}_{\mathrm{E}}{|}_i+{\alpha }_{10}{\mathbf{v}}_{10}{|}_i,$$

(7)

where the unitless parameter ${\alpha }_{10}$ is chosen in an ad hoc manner, typically in the range 0.01–0.03 (14), to measure windage, and ${\mathbf{v}}_{10}$ is the wind velocity at 10-m above the sea level. This is fundamentally different from the windage parameter α in the eBOMB model. Although the carrying flow velocity $\mathbf{u}$ in the eBOMB model, as described in Eq. 1, appears similar to the leeway model in Eq. 7, $\mathbf{u}$ is not arbitrarily chosen. Instead, it results from vertically averaging the Stokes drag exerted by the ocean velocity on the submerged part of an inertial particle and by the wind velocity on the exposed part (26). For example, we will refer to Eq. 7 with ${\alpha }_{10}=0.01$ as “Leeway 1%.” The accuracy of both the eBOMB and leeway models will be evaluated against Sargassum coverage distributions obtained from satellite data (see Materials and methods for more details on how these distributions are derived).

In this research, we focus our simulations on the year 2018, a time characterized by significant Sargassum presence in the Caribbean Sea. These simulations employ Eulerian ocean velocity (⁠${\mathbf{v}}_{\mathrm{E}}$⁠) derived from an amalgamation of various observational datasets, alongside Stokes drift (⁠${\mathbf{v}}_{\mathrm{S}}$⁠) and wind velocity (⁠${\mathbf{v}}_{\mathrm{W}}$⁠) generated by a reanalysis system. The parameters for our primary simulation were determined through a synthesis of heuristics, established values, and optimization. We initialize both the eBOMB and leeway models based on the Sargassum coverage distribution accumulated during the first week of March, as illustrated in Fig. 2, and subsequently simulate each model over a 14-week period. Details regarding the construction of the velocity fields, parameter optimization, and the computation of Sargassum coverage distribution are provided in Materials and methods.

Fig. 2.

AFAI-derived Sargassum coverage distribution from which the simulations in the subsequent figures were initialized. The value the colorbar is $-{\log }_{10}b$⁠, where b is the monthly fraction of Sargassum in a bin. Highlighted are regions encompassing Puerto Rico Island (solid rectangular box) and the Yucatan Peninsula (dashed rectangular box), chosen for model validation.

A qualitative analysis of the results for two representative weeks is illustrated in Fig. 3, where the eBOMB model is configured to use only physics. We consider Leeway 1% as it appears to be somewhat more accurate than Leeway 3% as per our more quantitative results shown in Fig. 4. The comparison indicates an overall enhancement achieved by the eBOMB model. Notice the wave-like pattern in the observed Sargassum distribution during the first week of April (Fig. 3, left column). This pattern is more accurately represented by eBOMB compared to Leeway 1%. Similarly, the Leeway 1% distribution near the Antilles Arc is shifted northwards relative to the observed distribution, while the eBOMB model captures its main location more precisely. The eBOMB model also shows better overall performance compared to Leeway 1% in the second week of May (Fig. 3, right column). Notably, Leeway 1% predicts high concentrations of Sargassum in the Gulf of Mexico, which the eBOMB model does not, aligning with observations. Additionally, Leeway 1% infers higher concentrations in the western Caribbean Sea and east of the Antilles Arc. The eBOMB distributions generally align more closely with observed distributions.

Fig. 3.

Comparison of the (C, D) eBOMB (with biology turned off) and (E, F) leeway models to (A, B) satellite-derived Sargassum coverage distributions for two different weeks in April 2018. The colorbar is as defined in Fig. 2, except that the scales in each week are different to maintain clarity of comparison.

Fig. 4.

Weekly JSD similarity between observed Sargassum coverage distribution and as predicted by the eBOMB model, with and without life cycles represented, and the leeway models, calculated in the region of the simulations.

A more rigorous quantitative comparison is achieved using the Jensen–Shannon Divergence (JSD) (38), a method for assessing the similarity between two probability distributions (see Materials and methods for details). The JSD is inherently non-negative, with a value of 0 signifying identical distributions, and values greater than 0 indicating divergence between the distributions. We commence by calculating the JSD between the observed Sargassum distribution and the predictions generated by the eBOMB model, both with and without biological factors, as well as Leeway 1% and Leeway 3%, across the entire simulation domain. This domain encompasses the Caribbean Sea, the Gulf of Mexico, and the North Atlantic from the equator to 32.5°N and west of 40°W. The results are presented adjusted for cloud coverage at weekly intervals in Fig. 4. Two observations are particularly noteworthy. First and foremost, the eBOMB model markedly outperforms both Leeway 1% and Leeway 3% (with the former being more accurate than the latter), providing robust validation for our proposed model. Secondly, the inclusion of biological aspects enhances the performance of the eBOMB model relative to observations, albeit this enhancement is marginal.

We further test the local accuracy of each model in two regions of interest, highlighted with boxes in Fig. 2. One region encompasses the Yucatan Peninsula, with its eastern coastline being notoriously affected by Sargassum chocking. The other region surrounds Puerto Rico Island, also reported to be affected by the accumulation of Sargassum, especially in its southern coasts. We again compute the JSD, renormalized in each of these regions by the total probability they contain. The results are shown on a weekly basis in Fig. 5. Note the overall better performance of the eBOMB model with no biology compared to Leeway 1% and Leeway 3% (Fig. 5, top panels). In bottom panels of Fig. 5, we restrict the comparison to eBOMB model with and without life cycles represented. As noted above, only a marginal overall improvement is provided by the activation of biology.

Fig. 5.

A, B) As in Fig. 4, but within two regions of interest, indicated with boxes in Fig. 2. C, D) As in (A, B), but restricted to the eBOMB model showing the effect of adding biology.

It is important to acknowledge that the quality of the results presented is contingent upon several factors. Regarding physical parameters, the representations of ocean and wind velocities are far from perfect, particularly the ocean velocity, which, despite being a data synthesis rather than a simulation outcome, omits critical aspects of ocean circulation such as submesoscale motions. Conversely, it has been observed that the Stokes drift component of ocean velocity exerts a negligible effect; however, there remains a possibility that such drift is inaccurately represented by reanalysis systems. Additionally, the wind velocity data provided at 10-m above sea level from reanalysis systems does not accurately represent the near-surface wind field. In any event, should an accurate operational coupled ocean–wave–atmosphere model become available, its output could be utilized to enhance the performance of the eBOMB model. The effects of wave breaking, for instance, may be parameterized in eBOMB by simulating Sargassum sinking in a manner analogous to the treatment of Sargassum mortality.

Concerning biological factors, it is postulated that biological effects may accumulate over time, thereby reducing discrepancies with observations. However, it is crucial to recognize that the nutrient and temperature data employed to inform the life cycles model Eq. 6 contain uncertainty. Moreover, our selection of the life cycles model is admittedly rudimentary. Enhancements are anticipated through the adoption of a more sophisticated model, which can be seamlessly integrated into the eBOMB framework.

Notwithstanding the aforementioned limitations, the eBOMB model, in its current configuration, surpasses the leeway model, which is further incapable of accurately depicting the observed aggregation of Sargassum within mesoscale eddies (32). Mesoscale eddies serve as significant mechanisms for the long-range transport of ocean properties (39–41), potentially contributing to the connectivity between distant regions of Sargassum concentration. Below, we demonstrate the capability of eBOMB particles to represent concentrations within mesoscale eddies.

Behavior near mesoscale eddies

In Ref. (27), a rigorous result pertaining to the dynamics of eBOM particles in the proximity of mesoscale eddies is presented. This result is equally applicable to immortal eBOMB particles, namely, those that are not influenced by physiological changes. An informal articulation of this result is presented below (see the Supplementary material for a rigorous statement and the corresponding proof).

 

Theorem 1.

In the presence of quiescent wind conditions, the cores of anticyclonic mesoscale eddies with coherent material boundaries function as finite-time local attractors for immortal eBOMB particles, whereas the cores of cyclonic eddies serve this role provided the particles exhibit sufficiently strong elastic connectivity.

By a coherent material boundary, we denote a loop formed by water parcels that experience an identical cumulative material rotation relative to the material rotation of the surrounding water mass (42). Numerical observations have demonstrated that each subset of these specialized material loops exhibits approximately uniform stretching over a finite-time interval (43). This attribute prevents them from experiencing breakaway filamentation (44). (See Supplementary material for a review of the notions of material coherence.) The property of immortal eBOMB particles articulated in Theorem 1 has significant implications for their long-range transport, and consequently for the transport of Sargassum rafts, as illustrated in the Caribbean Sea by (32). This property is not shared by leeway particles, which, under the conditions specified in Theorem 1, are unable to aggregate within mesoscale eddies, rendering these eddies ineffective as transport mechanisms.

The validity of Theorem 1 is evaluated in Fig. 6 (resp., 7) for the scenario of a mesoscale anticyclonic (resp., cyclonic) eddy. Both eddies were extracted from the same altimetry/wind/drifter synthesis of ocean velocity utilized in the aforementioned simulations. The extraction technique employed is geodesic eddy detection (44), ensuring that the extracted eddies possess material boundaries that exhibit uniform stretching over the extraction time interval. The red dots in each figure denote immortal eBOMB particles, while the blue dots represent leeway particles. The eBOMB particles in the top row of each panel have lower stiffness compared to those in the bottom row. Specifically, $A=15.1{\text{d}}^{-2}$ (resp., $A=200{\text{d}}^{-2}$⁠) in the top (resp., bottom) row. Further, the number of nearest neighbor connections has been increased, and the natural length of the spring has been recomputed as described in Materials and methods for this new initial distribution. The remaining eBOMB parameters are consistent with those used in the previous simulations, except that physiological transformations were considered. The leeway windage is set to ${\alpha }_{10}=0.01$⁠. Observe the tendency of immortal eBOMB particles to aggregate within the eddies. The expulsion of leeway particles is attributable to wind action, as observed during the simulation period (leeway particles with ${\alpha }_{10}=0$ represent water particles, which are unable to traverse coherent material boundaries; thus, if initially inside, they would remain enclosed). The outcomes of the tests align with the predictions of Theorem 1, despite being conducted beyond the strict domain of the theorem’s applicability, which is noteworthy.

Fig. 6.

Depictions of the temporal progression of immortal eBOMB particles (red dots) and leeway particles (blue dots) initially encompassing an anticyclonic mesoscale eddy, extracted on 2018 May 15 in the eastern Caribbean Sea from a synthesis of altimetry, wind, and drifter ocean velocity data utilizing geodesic detection with a coherence horizon of 2 months. The immortal eBOMB particles in (A, B, C) have lower stiffness compared to those in (D, E, F).

Fig. 7.

As in Fig. 6, but for the case of a cyclonic eddy.

Summary and outlook

In this study, we have developed a mechanistic model for pelagic Sargassum dynamics that integrates physical and biological processes. From a physical perspective, the model considers ocean currents and wind patterns to simulate Sargassum movement influenced by inertial effects. This was achieved by extending a Maxey–Riley theory for floating finite-size particles to incorporate nonlinear forces to simulate interactions within and between Sargassum rafts. Furthermore, the model includes a biological component by simulating Sargassum reproduction through vegetative growth and fragmentation, providing a comprehensive view of its life cycle. The model was validated using satellite-derived data, demonstrating superior performance compared to the leeway approach to transport, which is the foundation of current Sargassum modeling. Unlike our proposed model, the leeway model cannot simulate observed aggregation within mesoscale eddies, which provide a mechanism for long-range transport. Additionally, to encourage transparency and support broader scientific research, the model codes have been made publicly accessible, allowing the community to refine and apply the model in various research contexts, as well as in management and response efforts.

Future research endeavors will leverage the Maxey–Riley modeling framework to address critical questions regarding Sargassum connectivity among the Sargasso Sea, the tropical North Atlantic, and the Intra-America Sea. Furthermore, we will examine the potential for riverine runoff in the Gulf of Guinea to initiate blooms, which, under the influence of ocean currents and winds modulated by inertial effects, could result in elevated Sargassum concentrations throughout the tropical North Atlantic. This hypothesis necessitates validation through the application of suitable models and analytical tools. Notably, historical records indicate the presence of Sargassum in the Gulf of Guinea in 2009, predating the first significant invasion of the Caribbean Sea by Sargassum rafts (45), and preceding the extreme North Atlantic Oscillation (NAO) event hypothesized to have transported Sargassum from the Sargasso Sea to the tropical Atlantic (10). The model developed in this study is aptly suited for this purpose, and specialized methodologies from probabilistic dynamical systems (46) can facilitate the conceptualization of connectivity (47, 48).

Materials and methods

Sargassum distributions from AFAI

The Alternative Floating Algae Index (AFAI) is a metric derived from multiple spectral bands captured by the MODIS (Moderate-Resolution Imaging Spectroradiometer) instruments aboard NASA’s Aqua and Terra satellites, which can be utilized for the detection of floating algae, including Sargassum (49, 50). AFAI and its extensions have been used to create Sargassum coverage distributions with high precision and superior cloud rejection (51–53). Here, we begin with AFAI data accumulated in 7-day increments for the year 2018 (54, 55). We then apply a simplified version of data processing pipeline of (50) to produce weekly Sargassum coverage distributions, mainly taking the distribution to be given by the difference between local AFAI values and the median background in a certain pixel window, subject to thresholds.

Specifically, computation of a Sargassum coverage distribution from raw AFAI data proceeds in three steps which are shown in Fig. 8. First, the data are cleaned by removing values in the Pacific and data within 20 pixels of a coastal region. Then, each pixel is classified according to whether it contains Sargassum. To do this, the value of each pixel is decreased by the median value of all pixel’s surrounding it within a window of size 51 pixels. If the resulting value is greater than $1.79\times {10}^{-4}$⁠, it is considered to contain Sargassum. Finally, the distribution is created by scaling the AFAI values according to a linear interpolation between $-8.77\times {10}^{-4}$ and $4.41\times {10}^{-2}$ (recall that AFAI can take negative values), thresholding the top 85% of data to reduce noise and binning at roughly $1/2^{\circ }$ resolution. Bins with greater than half of the data missing are considered to be cloud-covered. The data are normalized to a probability distribution such that the total coverage for each month is equal to 1.

Fig. 8.

The main stages of the AFAI processing pipeline: A) raw AFAI data; B) Pacific and coastal data removed; C) locations of Sargassum-containing pixels in red; and D) the final thresholded and binned distribution, as in Fig. 2. Bins with significant missing data are colored in black.

Metric for model validation

The Jensen–Shannon Divergence (JSD) (38) is employed as the main metric to measure the difference between simulations and observations. Results obtained with other metrics, which exhibit similar patterns, are included in Supplementary material. Let $P_1$ and $P_2$ be two discrete probability distributions defined on a sample space $\mathcal{X}$⁠. Let $Q=\frac{1}{2}(P_1+P_2)$ be the mixture distribution of $P_1$ and $P_2$⁠, i.e. the cumulative distribution function of Q is equal to the sum of the cumulative distribution functions of $P_1$ and $P_2$⁠, each weighted by $\frac{1}{2}$⁠. Then, we have

$$\text{JSD}[P_1\Vert P_2]=\frac{1}{2}{\displaystyle \sum _{x\in \mathcal{X}}}{P}_1(x)\log \left(\frac{P_1(x)}{Q(x)}\right)+\frac{1}{2}{\displaystyle \sum _{x\in \mathcal{X}}}{P}_2(x)\log \left(\frac{P_2(x)}{Q(x)}\right).$$

(8)

We note that this is defined for one or both of the cases $P_1(x)=0$ and $P_2(x)=0$ by setting the term in the associated sum equal to zero.

Carrying flow representation

The characterization of the Eulerian ocean velocity ${\mathbf{v}}_{\mathrm{E}}$ is provided by an amalgamation of geostrophic flow, inferred from multisatellite altimetry data (56), and Ekman drift, induced by wind reanalysis (57). This synthesis is calibrated to align with the velocities of drifters from the NOAA Global Drifter Program (58), which are drogued at 15-m depth. The resultant altimetry/wind/drifter synthesis is available for distribution through ftp://ftp.aoml.noaa.gov/phod/pub/lumpkin/decomp. The product has been validated in experiments involving small object tracking in the Atlantic (28, 29). The Stokes drift velocity ${\mathbf{v}}_{\mathrm{S}}$ and the wind velocity at a 10-m elevation ${\mathbf{v}}_{10}$ are produced by the ECMWF Reanalysis v5 (ERA5), accessible via https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5. In our simulations we have taken ${\mathbf{v}}_{10}$ to be representation of the near-surface wind velocity ${\mathbf{v}}_{\mathrm{W}}$⁠.

Determination of parameters via optimization

In order to determine the free parameters of the eBOMB model, we optimize simulations by treating the AFAI-derived Sargassum distributions as a target. We initialize our clumps to match the Sargassum coverage distribution in the first week of March 2018. This is done by dividing the nonzero coverage values into 10 logarithmically spaced levels labeled $i=1,2,\dots ,10$ so that each longitude–latitude bin has an i associated to it by that level. Then, we place i clumps uniformly at random in each bin, resulting in 2,710 initial clumps. A simulation was then run for a time up to and including the second week of May 2018, recording the location of each clump every 0.1 days. Note that our time-series comparisons, e.g. Fig. 4 extend beyond this record. For the time between 2018 March 7 and 14, every clump position was pooled and binned to the same resolution as the AFAI-derived Sargassum distribution during this week. The simulated and satellite distributions were both normalized and the JSD was computed using Eq. 8, excluding bins that were covered by clouds. This process was repeated for each week in the full integration time span, and the total loss for the simulation is taken to be the sum of all weekly JSD values. In this way, we obtain a single quantity which, when small, indicates that the simulation is accurate relative to the observations in a cumulative sense. Then, for a given set of parameters to optimize, we use the cumulative JSD value as the target to minimize.

To obtain our final values in Table 2, we optimized in two stages. In the first stage, we turn off biological growth and death and optimize the physics parameters $\tau ,\delta ,\sigma$ and A. In the second stage, we fix the physics parameters, turn on biological growth and death, and optimize ${\mu }_{max},m,k_{\mathrm{N}},S_{min}$ and $S_{max}$⁠. This staged optimization ensures that the core physics model is accurate without relying too heavily on biological factors. We therefore have two derivative-free black-box optimization problems with relatively expensive function evaluations. To solve this, we use the NOMAD implementation of the mesh adaptive direct search algorithm (59), which interfaces with the SciML optimization packages. The NOMAD solver accepts box constraints; we set these according to our confidence in the value of each parameter. The buoyancy was measured in (60) to lie in the range $1.00$ to $1.49$⁠. The Maxey–Riley theory we apply requires the inertial response time to be small (see Supplementary material); it is sufficient to take $\tau f\lesssim 1$ to obtain an upper bound and the lower bound is a factor of 10 below this. Similarly, we take the Stokes drift coefficient in the range 0.5 to 1.5. Other optimized parameters were constrained within approximately a factor of 10 on either size of heuristic starting values. Each stage of the optimization was run with a hard time limit of 1 h on a 2022 MacBook Pro (Apple M2 chip, 16 GM RAM). We note that certain biological parameters in our study exhibit values considerably lower than those reported in the literature (14, e.g.). This discrepancy can be attributed to the incorporation of biological effects as a supplementary module to the mechanistic framework within our model, thereby introducing expected variations. Furthermore, our biological model is implemented in the context of discrete clumps, wherein the concept of fractional growth is less rigorously defined. We anticipate that the breadth and precision of the growth and death effects can be substantially enhanced in subsequent research, particularly to facilitate the direct inclusion of empirical measurements.

Complete model specification

Tables 1 and 2 summarize all of the datasets and parameters required to fully define our model and Algorithm eBOMB provides the pseudo-code. We comment briefly on an implementation detail regarding the life cycle portion of Algorithm eBOMB. The SciML ecosystem provides a callback interface where the integration status can be conditionally modified during or after each step in an efficient and safe manner. This allows us to log the clump positions at regular time steps $(t_i,t_i+h)$ while the underlying integration checks for growth and death (due to either biology or beaching) at the end of whichever step size the integrator has chosen. Hence, we can still take advantage of fast adaptive integration algorithms such as the method of (61). Further, using these callbacks, Algorithm eBOMB can be implemented such that the length of c, that is, the length of the vector of clump positions, is dynamically modified whereby at given time it only contains the positions of living clumps. This can provide an increase in speed at the expense of introducing another variable whose length is always equal to the length of c that keeps track of the labels of the clumps by vector position.

Table 1.

Model datasets.

Name	Symbol	Units	Source
Eulerian velocity	$v_{\mathrm{E}}$	${\text{km d}}^{-1}$	(62)
Stokes drift	$v_{\mathrm{S}}$	${\text{km d}}^{-1}$	(63)
Wind velocity	$v_{\mathrm{W}}$	${\text{km d}}^{-1}$	(63)
Temperature	T	°C	(64)
Nitrogen content	N	mmol m−3	(65)
Landmass location	land	None	(66)

Water and wind velocities are shared by eBOMB and Leeway.

Table 2.

Parameters of the eBOMB model.

Parameter name	Symbol	Value	Units	Source
Windage parameter	α	$3.37\times {10}^{-3}$	None	Derived
Inertial response time	τ	$0.0103$	d	Optimization
Buoyancy	δ	1.14	None	Optimization
Clump radius	a	$6.4\times {10}^{-5}$	km	Derived⋆
Spherical factor	R	0.823	None	Derived
Coriolis parameter	f	2.18213	d	Derived
Stokes drift coefficient	σ	1.2	None	Optimization
Stiffness amplitude	A	15.1	${\text{d}}^{-2}$	Optimization
Stiffness cutoff scale	Δ	0.2	km	Heuristic
Spring natural length	L	Variable	km	Heuristic
Maximum growth rate	${\mu }_{max}$	0.00541	${\text{d}}^{-1}$	Optimization
Mortality rate	m	0.00402	${\text{d}}^{-1}$	Optimization
Nitrogen uptake half saturation	$k_{\mathrm{N}}$	0.000129	mmol/m3	Optimization
Temperature growth minimum	$T_{min}$	10.0	°C	(12)
Temperature growth maximum	$T_{max}$	40.0	°C	(12)
Clump amount minimum	$S_{min}$	$-0.00482$	None	Optimization
Clump amount maximum	$S_{max}$	0.001	None	Optimization

The designation “Derived^⋆” signifies that a constitutes a more fundamental parameter in comparison to τ; nevertheless, τ is optimized due to its more explicit occurrence in the eBOMB model equations.

Algorithm 1

eBOMB model for Sargassum clump trajectories with life cycles.

Require: Datasets from Table 1 and constants from Table 2.

Require: Integration time span $(t_1,t_2)$⁠.

Require: Integration algorithm with time step h. ⊳ Such as RK4.

Require: The maximum number of allowed clumps $N_{max}$⁠.

Require: Clump initial positions $c=(x_1,y_1,x_2,y_2,\dots ,x_{N_{max}},y_{N_{max}})$⁠.

Require: Clump initial connections (neighbor(i, ${t_1))}_{1\le i\le N_{max}}$⁠.

Require: Clump initial amounts $S=(0{)}_{1\le i\le N_{max}}$⁠. ⊳ Or a distribution of S values.

Require: Clump initial life status $\ell =({\ell }_1,{\ell }_2,{\ell }_3,\dots {\ell }_{N_{max}})$ where ${\ell }_i\in \{1,0\}$ according to whether clump i is alive.

1: Compute L using c and the procedure described in Model overview with $K=5$⁠.

2: Set $t\leftarrow t_1$⁠.

3: Set $N\leftarrow \text{sum}(\ell )$⁠. ⊳ This is the total number of clumps that have ever existed.

4: While  $t<t_2$  do

5:  Integrate $\dot{c}=g(t,c)$ with g implied by Eq. 2 for one time step h.

6:  Set $c\leftarrow c(t+h)$ and $t\leftarrow t+h$⁠.

7:  Set ${\ell }_i\leftarrow 0$ for all clumps for which $(x_i,y_i)\in \text{land}$⁠.

8:  for  $\{i:{\ell }_i=1\}$  do

9:   Update $S_i$ using Eq. 5 evaluated at $(t,x_i,y_i)$⁠.

10:   If  $S_i<S_{min}$  then

11:    Set ${\ell }_i\leftarrow 0$⁠.

12:   else if  $S_i>S_{max}$  then

13:    Set $S_i\leftarrow 0$⁠. ⊳ Refresh this “parent” clump.

14:    Set $N\leftarrow N+1$⁠.

15:    Set ${\ell }_N\leftarrow 1$⁠.

16:    Generate $\theta \in [0,2\pi )$ uniformly at random.

17:    Set $x_N\leftarrow x_i+L\text{cos}\theta$ and $y_N\leftarrow y_i+L\sin \theta$⁠.

18:   end if

19:  end for

20:  Form new connections (neighbor(i, t) : ${\ell }_i=1$⁠) among living clumps.

21: end while

22: Return c

Acknowledgments

We extend our gratitude to Nathan Putman for the benefit of many useful discussions regarding Sargassum phenomenology.

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

This research was funded by the National Science Foundation (NSF) under grant number OCE2148499.

Author Contributions

G. Bonner: Conceptualization; Formal analysis; Code development; Writing—main text. F. J. Beron-Vera: Conceptualization; Writing—supplementary material, review & editing; Funding acquisition. M. J. Olascoaga: Conceptualization; Writing—review & editing; Funding acquisition.

Preprints

A preprint of this article is available at https://doi.org/10.48550/arXiv.2410.01468.

Data Availability

All data and initial source code are available through https://doi.org/10.6084/m9.figshare.27146685.v2 and further code developments will be available at the GitHub repository https://github.com/70Gage70/Sargassum.jl. All code utilized for the simulation, optimization, and visualization of the eBOMB model is disseminated freely in an open-source manner under the MIT License via the GitHub repository https://github.com/70Gage70/Sargassum.jl. The main tool for end users is the Julia package Sargassum.jl, distributed initially as version 0.1.6. This is a meta-package collecting all of the required sub-packages in one place, including functionality for simulations, interfacing, and AFAI data processing. In particular, we provide an interface that allows our software to be used with zero coding in the form of a reactive Pluto.jl notebook (67). All packages were developed using Julia 1.10.2 (68). All raw data corresponding to each field in Table 1 employed in this study are readily accessible for download through the provided interfaces. The integration of ordinary differential equations and black-box optimization interfaces is facilitated by the SciML ecosystem (69). The data visualizations presented in this article were generated using the Makie ecosystem (70).

Author notes

Competing Interest: The authors declare no competing interests.