Authors: Xing Zhou, Lyuba Novi, Mark E. Hay, Joseph P. Montoya, Aderinsola Aliu, Matthew J. Realff, Annalisa Bracco
Categories: Article, Marine biology, Biooceanography
Source: Nature Communications
Authors: Xing Zhou, Lyuba Novi, Mark E. Hay, Joseph P. Montoya, Aderinsola Aliu, Matthew J. Realff, Annalisa Bracco
Inundations of pelagic Sargassum plague the tropical Atlantic, with size and impacts steadily increasing to surpass 30 million tons in 2025. Understanding the drivers of Sargassum growth in the so-called Great Atlantic Sargassum Belt is fundamental to developing effective mitigation strategies for affected nations. We present a nonlinear regression model that both explains the seasonal and interannual variability observed between 2011 and 2022 and predicts Sargassum concentrations in 2023 and 2024. The growth of Sargassum, initiated by a prolonged negative phase of the North Atlantic Oscillation, is initially enhanced through winter mixed layer deepening in response to stronger winds. An additional overlooked driver is the recycling of nutrients within the mixed layer, carried out by the community of organisms associated with Sargassum and aging Sargassum mats. This contribution increases over time to become dominant in recent years, offsetting the increase in stratification in 2023 and 2024.
Massive and increasing inundations of pelagic Sargassum spp. (mostly S. natans and S. fluitans, referred to simply as Sargassum herein) have plagued the shores of the Caribbean and Gulf of Mexico since 2011^1^. Historically found in the Sargasso Sea and endemic to the subtropical Atlantic Ocean, Sargassum blooms have been developing in recent years also across the Intra-Americas Sea and the tropical North Atlantic, from West Africa to South America, where its growth is usually bounded by the South Equatorial Current and the North Equatorial Counter Current (Fig. 1). This so-called Great Atlantic Sargassum Belt or GASB has surpassed 20 million tonnes of Sargassum and stretched for over 8000 km at its monthly peak, usually June or July, every year since 2018. Given its large-scale coverage, growing trend, and associated major remediation costs, the GASB represents an expensive and ecologically worrisome challenge for Caribbean and other nations with limited resources.Fig. 1The GASB.The Great Atlantic Sargassum Belt (GASB) in July 2023. NASA Scientific Visualization Studio (https://svs.gsfc.nasa.gov/5298/).
The GASB was initiated by a prolonged, highly negative phase of the North Atlantic Oscillation during 2009–2010, which caused anomalous southward winds that transported Sargassum spp. from the Sargasso Sea to about 5 °N in the eastern Atlantic^2,3^. Its yearly recurrence and especially its intensification, however, remain a conundrum, begging the question of nutrient supply to support Sargassum growth. In the absence of a robust conceptual understanding of the drivers of the GASB growth, the capacity to predict its long-term evolution and to formulate cost-effective strategies for mitigating Sargassum inundations remains constrained.
The fast growth of the GASB in recent years has been variously attributed to increased nitrogen runoff from the Amazon and Congo Rivers, increased atmospheric deposition, increased coastal upwelling, increased vertical ocean mixing, increased sea surface temperatures (SST), and an equatorial upwelling of phosphorus driving nitrogen fixation^1–8^. While there has been an increase in nitrogen discharge by the Amazon between 2014 and 2018, a direct, annual correspondence between riverine NO3^−^ flux and Sargassum biomass is not apparent in the Hydrology and Geochemistry of the Amazon basin (HYBAM) observatory data^8,9^. Coastal upwelling cannot explain the offshore genesis of the GASB, while atmospheric deposition at the required rate would be detectable. More recently Podlejski and co-authors^6^ pointed to surface ocean temperatures to explain growth and decay of the blooms, but although temperature may explain variations in the annual cycling of recent blooms, temperature alone cannot support the extraordinary growth trend, nor the nitrogen and phosphorus enrichment of the GASB Sargassum populations compared to those in the Sargasso Sea habitat^7^. Lastly, equatorial upwelling^8^ has declined since 2022, while Sargassum concentrations have continued to increase.
Among the proposed mitigation approaches, much interest revolves around exploiting the Sargassum blooms as a marine carbon dioxide removal (mCDR) strategy by either sinking the macroalga or converting it into biofuels by leveraging the fact that carbon comprises ~27–30% of the dry mass of Sargassum collected in the intra-America Seas^9^. Conversion can generate solid, liquid, or gaseous fuel. Solid fuel methods include direct combustion, which is simple but produces poor-quality fuel and harmful emissions, torrefaction, which heats Sargassum in the absence of oxygen to produce charcoal with higher energy density, and densification, where Sargassum is compressed into bio-pellets for efficient fuel use in stoves and boilers^10^. Liquid fuel processes (i) fermentation, where sugars from Sargassum are converted to ethanol by the yeast Saccharomyces cerevisiae; (ii) pyrolysis, which produces bio-oil and gases at high temperatures; (iii) liquefaction which involves the use of biomass at low temperature and high pressure to form bio-oil; and (iv) transesterification, where lipids from the Sargassum are extracted and converted into biodiesel^11–13^. In gaseous fuels, gasification converts Sargassum into producer gas (CO, H₂, CH₄) for energy use, while anaerobic digestion produces biogas (methane) through microbial decomposition^14^. Each method varies in energy efficiency, environmental impact, and practical use but in all cases the scalability of mCDR options depends on the drivers of the GASB and on its likelihood of continuing to plague the tropical and subtropical North Atlantic.
In this work, we investigate how nutrients are supplied to the GASB, why this supply has increased in the past decade, and we build a predictive model of past and future Sargassum concentrations. Our nonlinear regression model describes the growth of the GASB from January 2011 to December 2022, considering observational data and an ecological hypothesis. Its robustness is tested by predicting the observed Sargassum concentrations of 2023 and 2024. Specifically, we show that the evolution of the GASB can be modeled and the Sargassum concentrations predicted with a high degree of accuracy as a response to mixed-layer deepening, and consequent upwelling of nutrients, with an additional nutrient source provided by the community of organisms associated with the Sargassum. The contribution of mixed-layer deepening has been decreasing since 2018, with the input from the community of organisms within and around Sargassum consistently growing over time to dominate the nutrient supply to the GASB in the past five years.
We built a non-linear regression model to describe the monthly evolution of Sargassum biomass concentration from 2011 to 2022. We then used the model to predict concentrations in 2023 and 2024.
First, we selected the variables that showed a statistically significant change (see “Methods”) over the latitudinal band 1–15°N in the direction of influencing and potentially amplifying Sargassum growth in the GASB. We tested mixed-layer depth (MLD), SST, sea surface salinity (SSS), eddy kinetic energy, eolian dust deposition, and performed a hotspot of change evaluation^15^, in each season separately and for yearly data. This methodology allows consideration of changes in mean values, multi-year seasonal variability, and extremes. A correlation analysis was adopted for the time series of major riverine flow into the Tropical Atlantic. Among the variables, only MLD and dust satisfied our criterion of statistical significance (Fig. 2 and Supplementary Fig. 1). Mixed-layer deepening would allow for greater entrainment of nitrogen and phosphorus-rich mode waters into the surface mixed layer^16^, while an increase in dust from the African continent would supplement micronutrients. SST was also retained despite showing only a modest correlation and exerting limited influence on model skill, as warming may improve the conditions for Sargassum growth^5^. Riverine input of nitrate (NO3^−^) did not show any significant correlation with Sargassum biomass, in agreement with^8^ (Supplementary Fig. 2).Fig. 2Mixed-layer deepening in the tropical Atlantic.A1–A4 (Left): Maps of mixed layer depth change expressed as percentage differences between the periods 2011–2022 and 1999–2010, winter to fall (top to bottom). A 20% change inside the GASB corresponds to roughly 8 m in MLD. B1–B4 (Right): Time series of MLD anomalies over 1999–2022 over the GASB area (black line in left panels) for each season. Source data are provided as a Source data file.
Our analysis revealed that the mixed layer deepened in the 2011–2022 period over the GASB area compared to the 1999–2010 time interval, despite an increase in near-surface stratification controlled by SST and salinity^17,18^. MLD changes were especially relevant in winter and spring and from 2012 to 2020, when part of the GASB domain deepened by more than 20%. An increase in Dust Aerosol Optical Depth (DAOD) input characterized the same season. SST changes, on the other hand, warmed the GASB by less than 0.2 °C with a nearly homogenous signal in winter and fall (Supplementary Fig. 3).
The simultaneous MLD deepening and dust aerosol optical depth (DAOD) increase point to wind intensification as the primary cause, as verified by exploring changes in wind intensity and wind stress (Supplementary Fig. 4). This trend is, in part, a response to interannual variability in the Atlantic Meridional Oscillation (AMO)^8^, and in part the result of a global intensification of wind energy in the second decade of the XXI century, in agreement with a recent analysis of land surface wind speed^19^. The intensification appears especially strong in the Canary upwelling region with more nutrient-rich waters being pushed towards the GASB from the northeast in winter, spring, and fall. This change in the wind field in an upwelling system follows a general trend in coastal upwelling systems associated with an increase in land-sea temperature contrast^20^.
In addition to terms accounting for MLD deepening, dust increase, and SST changes, our regression model included a biological contribution to describe the self-fertilization potential of the GASB. This term is proportional to the pre-existing biomass concentration under the hypothesis that old Sargassum mats and the Sargassumsphere - the community of interacting organisms associated with Sargassum that includes a diverse array of aquatic species- may contribute to the collection, concentration, and recycling of nutrients within the ocean mixed-layer, effectively harvesting nutrients from a broader volume around the Sargassum mat and concentrating nutrients to feed the bloom. Just as birds roosting on islands forage broadly for food, concentrate those nutritious resources back on the island, and by doing so alter the growth, composition, and nutrient levels in island vegetation and coastal oceans^21–23^, plankton feeding fishes, shrimps, crabs, and other invertebrates associated with Sargassum mats are thought to harvest plankton from surrounding waters and concentrate these nutrients in the Sargassum mats among which they shelter and excrete nutrients^9^. Fishes associated with Sargassum floats release nutrients^24^, and Sargassum can efficiently take-up short nutrient pulses and grow over considerable periods following the uptake^25^. When Sargassum decomposes, the remineralized nutrients become available for phytoplankton uptake. The remaining phytoplankton biomass and dissolved nutrients in the mixed layer continue to sustain the planktonic community as well as the Sargassum mats. These animals feed on passing phytoplankton and, through continual excretion, release nutrients back into the Sargassumsphere, effectively preconditioning the environment for the next bloom. With the surface Tropical Atlantic and specifically the GASB being among the most stratified regions due to salinity (e.g., ref. ^26^), the loss of dissolved nutrients out of the mixed layer through vertical mixing is slow, further contributing to their retention.
To support the Sargassumsphere hypothesis, Fig. 3 shows the nitrogen isotopic composition (δ^15^N) of Sargassum and of the common mobile epibionts, Latreutes fucorum, mixed L. fucorum and Leander tenuicornis (“shrimp”), and the crab Portunus sayi, all found in abundance in Sargassum mats. Samples were collected in coastal waters of the US Virgin Island of St. Thomas in May 2024. The δ^15^N of our Sargassum samples is lower than the δ^15^N of two major potential sources of nitrogen in our study region. Subsurface NO3^−^ typically has a δ^15^N of ca. 4.5‰ (e.g., refs. ^27–29^) while the δ^15^N of diazotrophs common in the tropical and subtropical Atlantic (Trichodesmium and Diatom-Diazotroph Associations, DDAs) typically ranges between −2 and 0‰^30–32^. In contrast, the δ^15^N of NH4^+^ excreted by animals is ca. 3‰ lower than the δ^15^N of the organic matter being catabolized^33,34^, or ca. 6‰ lower than the animal biomass. The common mobile epibionts on Sargassum have δ^15^N values ranging between 1 and 5‰, which implies an excretory flux of NH4^+^ with a δ^15^N of −5 to −1‰, low enough to account for the Sargassum δ^15^N values we measured.Fig. 3Nitrogen composition in support of the Sargassumsphere hypothesis.A Nitrogen isotopic composition (δ ^15^N) of Sargassum (sample number n = 10) and the common mobile epibionts Latreutes fucorum (n = 24), mixed L. fucorum and Leander tenuicornis (“shrimp”, n = 63), and Portunus sayi (n = 29). B Schematic illustrating the importance of excreted ammonium in producing the δ ^15^N of Sargassum, which is lower than the δ ^15^N of both subsurface nitrate and nitrogen recently fixed by diazotrophs. Source data are provided as a Source data file. A has been generated using R v4.5.2. B has been generated using Procreate v5.4, Adobe Illustrator v29.6.1, and Microsoft Powerpoint v16.101.
For this self-fertilization term, we hypothesize that nutrients are released in the mixed layer when the Sargassum bloom declines in fall and winter, and that in any given month, the nutrients available through this mechanism depend broadly on the 3-month average biomass concentration of the previous year (assumed to be null prior to 2011). We assumed, for example, that the Sargassum concentration in June 2018 depends on the average concentration observed over May-June-July (MJJ) of 2017. Given that the lifespan of the species abundant in the Sargassumsphere varies between 5 and 6 months (Latreutes fucorum) to few years (Portunus sayi), the annual timescale is sensible. Our choice of using the 3-month average concentration of the previous year can be rectified to instead adopt the annual average concentrations of the previous year without statistically significant changes in skill.
The Sargassumsphere hypothesis also assumes that not all Sargassum is lost to sinking in the GASB. According to^35^ Sargassum loss through sinking will occur only when the depth that a water parcel originally at the surface approaches or exceeds 100 m, which, in the tropical Atlantic, corresponds to strong mixing events occurring mostly in winter. This would support mats spending long times near the surface especially between spring and fall^36^. The Sargassumsphere contribution is further modulated by an attenuation or amplification factor linked to changes in mixed layer depth. For example, if the MLD in MJJ 2017 was deeper than average for that season, the concentration of nutrients remineralized in the mixed-layer because of the Sargassumsphere would be more diluted and vice versa.
In summary, the modeled Sargassum concentration is described by Eq. 2 (see “Methods”). The correlation coefficient (R) and mean square error (MSE) between the observed time series and the modeled ones are R = 0.84 and MSE = 4.34 (Fig. 4A).Fig. 4The regression model, key contributors to Sargassum growth, and its prediction skill.A Observed (black) and modeled (green) time-series of Sargassum concentrations in million of tonnes from Jan 2011 to Dec 2022. B Time evolution of two major factors influencing Sargassum mixed-layer deepening and self-fertilization. C Prediction of Sargassum concentrations obtained using Eq. 2, excluding the dust and SST terms. The model was trained on data from 2011 to 2022 and MLD knowledge is retained up to 3 months prior. Correlation coefficient R and mean square error MSE are calculated over the entire 2011–2024 period. For (A, C): DoF (degrees of freedom) = 43 considering the Sargassum data autocorrelated over 3 months and 5 parameters (144/3 – 5 = 43); p < 0.01. Source data are provided as a Source data file.
Sensitivity tests of the non-regression model were performed eliminating each term in Eq. (2). The dominant terms were found to be mixed-layer deepening and self-fertilization (Supplementary Fig. 5), and their relative evolution through the years is shown in Fig. 4B. This plot clearly shows that mixed-layer deepening was crucial to the blooms between 2011 and 2018, while the self-fertilization term is essential to capture the massive amounts of Sargassum recorded after 2018.
If MLD deepening has been key in sustaining the GASB, then its contribution would be reflected in the nutrient supply. We therefore calculated a time-series of nitrogen (N) input based on the MLD changes accounting for the lag that maximizes the correlation between Sargassum concentration and MLD (lag1 = −3 months). We used a reanalysis (see “Methods”) to evaluate how much excess N may have been available in the mixed layer over the GASB region in each month between January 2011 and December 2022, compared to the mean value of N for the same month over the 1999–2010 period (Supplementary Fig. 6). Considering that the N content is about 1.2% of the tissue elemental composition in GASB samples^9^, the amount of nitrogen introduced by MLD deepening is approximately 50% of that required to sustain the observed biomass in the 12 years considered.
Lastly, the predictive skill of our model was evaluated by testing its capacity to reproduce the Sargassum concentrations observed in 2023 and 2024. We used the coefficients derived from the 2011 to 2022 period alone, setting to zero those related to the dust and SST terms for the predicted years, and verified that the prediction skill of Sargassum inundations remains high whenever previous year concentrations and MLD three months prior are known (Fig. 4C). This is the case for years of extreme warming in the Atlantic, during which the seasonal cycle was been partially impacted, the MLD reverted to climatological values, and the Sargassum growth occurred earlier^37^.
If we can predict the amount of Sargassum that will accumulate over the next three months, we can proactively assess whether such an increase will lead to beaching using available forecast models for ocean currents^3,38^, and it would be possible to plan for offshore harvesting to prevent Sargassum from reaching the shore. This planning would involve determining the necessary number of laborers, and associated operational costs required for efficient harvesting. Once harvested, the Sargassum can be valorized into biofuels or other products, converting the burden into a valuable resource. Moreover, with accurate predictions, we can estimate the potential biofuel production per month depending on the available technologies. Our results show, however, that a sustainable harvesting process is critical to allow the natural replenishment of Sargassum and ensure long-term balance in the ecosystem.
Nutrient availability is key to the development of Sargassum populations^39^, and enhanced nutrient input is required to explain the sustained growth of the GASB since 2011^2^. This work shows that Sargassum growth in the GASB, with its strong interannual variability, was initially enhanced through mixed layer deepening with a significant, and increasing through time, contribution provided by the recycling and remineralization of nutrients in the mixed layer by the community of organisms associated with Sargassum and old Sargassum mats. Consequently, the recent increase in stratification of the tropical Atlantic in 2023 and 2024^37^, has not limited the blooms, since the Sargassumsphere contribution is now dominant. The schematic in Fig. 5 summarizes the drivers of primary productivity before 2011, between 2011 and 2020, and from 2020 onward in the Tropical Atlantic.Fig. 5Conceptual diagram of Sargassum growth.Drivers of primary productivity including uptake of new nitrogen (green arrows), trophic and export fluxes (black arrows), remineralization fluxes (red arrows), and nutrient injection into the mixed layer of the GASB region before 2011, between 2011 and 2020, and after 2020. The diagram has been generated using Procreate v5.4, Adobe Illustrator v29.6.1, and Microsoft Powerpoint v16.101.
The Sargassumsphere hypothesis and its role in self-fertilization of the blooms is supported by our isotopic measurements. The common mobile epibionts found in the Sargassum have δ^15^N values ranging between 1 and 5‰, which implies an excretory flux of NH4^+^ with a δ^15^N of −5 to −1‰, low enough to account for the Sargassum δ^15^N values we measured (Fig. 3).
Our nonlinear regression model represents a fundamental step towards understanding how the Sargassum blooms that have plagued the tropical Atlantic since 2011 are not just maintained but have been growing over time. Our model quantifies the seasonal predictability of these algal inundations and improves their seasonal forecast. In doing so, it offers societal value by opening new avenues for proactive planning and response, and by helping frame management responses to mitigate the problem.
The role of mixed layer deepening, especially relevant in winter and associated with increasing wind and wind stress intensity, suggests that climate variability drove the GASB amplification since its inception in 2011 up to about 2020. Mixed layer deepening has supplied about half of the nitrogen required by the GASB, with biomass from, or associated with, previous blooms contributing further nutrients above the base of the mixed layer, especially since 2020. The relative dominance of these two terms in explaining the observed tendencies ensures that the Sargassum amplification in the tropical Atlantic is likely to continue, and that the life cycle assessment of mCDR options can assume at least as large if not larger concentrations of Sargassum in the near future.
Monthly MLD, SSS, SST, and surface geostrophic velocity data from January 1999 to December 2022 were obtained from the CMEMS global ocean eddy-resolving GLORYS12V1 reanalysis at 1/12° horizontal resolution (downloaded using E.U. Copernicus Marine Service Information, 10.48670/moi-00021, accessed on 11-03-2026). We used the multi-year (MY) product available until June 2021 and the interim multi-year (MY-INT) product from July 2021 to December 2024. For MLD, the monthly ocean mixed layer thickness defined by sigma theta (mlotst variable) was considered. A comparison with MLD data from state-of-the-art reanalysis datasets at 0.25° horizontal resolution, the ORAS5^40^ (downloaded using E.U. Copernicus Marine Service Information, accessed on 11-03-2026) and SODA Version 3.3.1^41^ (downloaded from APDRC at https://apdrc.soest.hawaii.edu/datadoc/soda_3.3.1.php, accessed on 11-03-2026), was also performed.
The DAOD data were retrieved from the CAMS global reanalysis (EAC4)^42^ monthly averaged fields from January 2003 (first available month) to December 2022 at horizontal resolution of 0.75° (duaod550 – single level variable) (downloaded from the Copernicus Atmosphere Monitoring Service (CAMS) Atmosphere Data Store, 10.24381/d58bbf47, accessed on 11-03-2026).
For the calculation of multi-year changes the variables were bilinearly remapped to 1/4° horizontal resolution using CDO (Climate Data Operator) available at https://code.mpimet.mpg.de/projects/cdo, to isolate variability at the mesoscale and larger.
In addition, we retrieved the horizontal (U) and meridional (V) wind components at the 1000 hPa pressure level from the ERA5 monthly averaged reanalysis^43^ from Jan 1993 to December 2022, available at 0.25° horizontal resolution (downloaded from the Copernicus repository, Copernicus Climate Change Service (C3S) Climate Data Store (CDS). 10.24381/cds.adbb2d47, accessed on 11-03-2026).
For the evaluation of the nitrogen input associated to MLD deepening, we used the Global Ocean Biogeochemistry Hindcast (version GLOBAL_MULTIYEAR_BGC_001_029, downloaded using E.U. Copernicus Marine Service Information, 10.48670/moi-00019, accessed on 11-03-2026) based on the PISCES biogeochemical model and forced by the FREEGLORYS2V4 ocean physics, and the variable monthly mole concentration of nitrate in sea water (NO3) in [mmol/m^3^].
For the isotopic analysis shown in Fig. 3, we collected samples of Sargassum and common mobile epibionts from waters south of St. Thomas, USVI, in May 2024. All samples were frozen (−20 °C) immediately after collection and dried after transport to Atlanta. We measured the isotopic composition of our dried samples by continuous-flow isotope ratio mass spectrometry using a Micromass 100 interfaced to a Carlo Erba NA2500 elemental analyzer for online combustion and purification of sample nitrogen and carbon. We used both elemental (methionine) and isotopic (peptone) standards to check instrument stability and to correct for analytical blanks^44^. We conservatively estimate that the overall analytical precision of our isotopic measurements is better than ±0.1‰.
To select the variables included in the regression model, we first investigated which changed in a statistically significant way after 2011 compared to the earlier period through a hotspot of change analysis^15^. The non-parametric method, which does not assume a specific probability distribution for the data, is flexible and can be applied to datasets regardless of their distribution. This analysis consists in calculating a Standard Euclidean Distance index (SED) that aggregates the changes in means, variability and extremes of the variable being examined point-by-point according 1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {SED}=\sqrt{{\sum }{i=1}^{N\Delta }{\sum }{j=1}^{4}{\left(\frac{{\Delta }{{ij}}}{p95\left(|{\Delta }{{ij}}|\right)}\right)}^{2}}
Here *N*Δ is the total number of indicators per each variable, *i* the index identifying each indicator, *j* identifies the season, and p95 is the 95th percentile computed spatially considering all grid points. Therefore Δ~*ij*~ is the *i*th indicator in the *j*th season (December-January-February and so on). For each variable we considered (MLD, SST, SSS, dust, eddy kinetic energy) we evaluated changes in means, variability and extremes between two periods of equal length 2011–2022 and 1999–2010. For MLD and dust changes were statistically significant in both means and extremes at least in some season, for SST only in mean, while no significant changes were found for salinity and eddy kinetic energy. ### Multi-year seasonal means changes calculation To evaluate changes, we first compared the multi-year seasonal means over two periods of equal length, after the *Sargassum* bloom initiation in 2011 (2011–2022) and before it (1999–2010), separately in each season. For example, winter changes in MLD (Δ~MLD, DJF~) are computed as Δ~MLD, DJF~ = *ysm*(MLD~DJF~)~(2011-2022)~ − ysm(MLD~DJF~)~(1999–2010)~, where ysm is the multi-year seasonal mean of MLD and then expressed as percentual changes with respect to the pre-2011 conditions as Δ%~MLD, DJF~ = 100 × (Δ~MLD, DJF~/ysm(MLD~DJF~)~(1999–2010)~). ### Regression model For the regression model, we calculated the monthly time-series of MLD, DAOD, and SST after removing the seasonal cycle over the study area grid point by grid point. The de-seasonalized time-series were then spatially averaged over an area broader than the GASB region and bounded approximately between [89°W–15°W] and [1°N–15°N], to account for the horizontal nutrient transport from the surrounding ocean into the GASB domain. Other areas were tested as well (see Sensitivity section further below). For MLD and SST a running mean (5-months for MLD and 3-months for SST) was applied to the timeseries to remove high frequency variability, while for DAOD a 3-months moving sum accounted for its sliding cumulative effect (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_{{ms}} $$\end{document}Dms). We then determined the time lag within a range [−12 months, 6 months] that maximized the absolute value of Pearson’s linear correlation coefficient between the de-seasonalized time series of the anomalies (with respect to their mean value) of these averaged quantities and that of monthly *Sargassum* biomass (S~obs~). The lag is non-zero only for the MLD, for which *lag1* = −3 months with <MLD> preceding S~obs~ where < > indicates spatial averaging. The model is robust to the choice of a running mean between 2 and 5 months, and of a lag between -1 and -4 months. The functional choice of the MLD term accounts for the nonlinear increase of nutrient concentrations with depth. A simple linear dependence was chosen for SST, which anyway showed small changes, and dust, where it accounts for its accumulation. We modeled the self-fertilization term that depends on prior concentration of *Sargassum*, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {C}_{{lag}} $$\end{document}Clag, considering a cubic dependence on the MLD. Its functional form is therefore proportional to previous concentrations and inversely proportional to the volume where the concentrations retained close to the surface are easily diluted (i.e., the mixed layer). The modeled *Sargassum* biomass concentration, Y, is therefore described by the equation2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Y={fH}[{{\mathrm{sgn}}}\left(f\right)],{with\; f}=[{{{\mathrm{sgn}}}\left({MLD}\right)b}_{1}{|ML}{{D|}}_{{rm}-{lag}1}^{\left|{b}_{2}\right|} \\+{b}_{3}\left(1-{{MLD}}_{{rm}-{lag}2}^{3}\right){C}_{{lag}}+{{b}_{4}D}_{{ms}}+{b}_{5}{{SST}}_{{rm}}] $$\end{document}Y=fH[sgnf],withf=[sgnMLDb1∣MLD∣rm−lag1b2+b31−MLDrm−lag23Clag+b4Dms+b5SSTrm]where H(.) is the Heaviside function and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {{\mathrm{sgn}}}(.) $$\end{document}sgn(.) the sign function. The coefficients b~1~-b~5~ are determined by fitting Y to *S*~*obs*~. We also explored the use of a machine learning approach that chose the terms from a library of functions in an unsupervised manner. This approach, however, did not produce an easily explainable solution and was abandoned. ### Relative role of the regression model’s terms We quantified the relative role of each term in Eq. 2 by removing each contribution, one at a time, as shown in Supplementary Fig. 5. Removing the terms that depend on the MLD caused the correlation between model and observed time-series to decrease to *R* = 0.66. An even greater drop in skill was found when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {C}_{{lag}} $$\end{document}Clag (the Sargassumsphere contribution) was not considered (*R* = 0.58). On the other hand, only minor changes are induced by neglecting the dust \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ ({{b}_{4}D}_{{ms}}) $$\end{document}(b4Dms) or SST \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ ({b}_{5}{{SST}}_{{rm}}) $$\end{document}(b5SSTrm) terms, with *R* = 0.81 or *R* = 0.84, respectively. ### Sensitivity of the regression model #### Parameter sensitivity There are two parameters that we initially chose arbitrarily and that may have influenced the robustness of the non-linear regression model results. These are the moving window applied to the running mean for MLD (5 months) and SST (3 months), and moving sum for dust (3 months), and the time lag parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {lag}1 $$\end{document}lag1 for MLD. We defined \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {lag}1 $$\end{document}lag1 within a range from −12 months to +6 months, selecting the value that maximized the absolute value of the Pearson’s correlation coefficient between the averaged environmental variables and monthly *Sargassum* biomass. We found that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {lag}1 $$\end{document}lag1 values were consistently high between −1 and −4 months, with −3 months typically showing the strongest correlation. To evaluate the sensitivity of the model to these choices, we conducted sensitivity tests by adjusting the moving window from 2 to 5 months and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {lag}1 $$\end{document}lag1 value from −1 to −4 months. We compared the new model outputs to the model results presented in the main text (green area in Fig. 3) using two statistical metrics, the correlation coefficient (R) and the mean squared error (MSE). We tested the sensitivity separately by (1) changing the moving mean window for SST, (2) changing the moving sum window for dust, and (3) changing both the moving mean window and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {lag}1 $$\end{document}lag1 for MLD. All statistical metrics indicate that our non-linear regression model is extremely robust. For the SST sensitivity tests, in all cases *R* > 0.99 and MSE < 0.01. For dust, we obtained *R* > 0.98 and MSE < 0.40, and for MLD, *R* > 0.90 and MSE < 2.0. Furthermore, the decrease in R values and the increase in MSE values from SST to MLD suggest that *Sargassum* blooms are more sensitive to variations in MLD, which is consistent with the conclusions presented in the main text. For the Sargassumsphere term, its functional form should be proportional to previous *Sargassum* concentrations and inversely proportional to the volume where the concentrations retained close to the surface are diluted (i.e., the mixed layer). If so, from a mathematical standpoint, we can write an optimization problem in the 3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {b}_{3}{\sum }_{{lag}=0}^{{lag}=\tau }{W}_{{lag}}\times {C}_{{lag}}(1-{{MLD}}_{{lag}}^{3}) $$\end{document}b3∑lag=0lag=τWlag×Clag(1−MLDlag3)where each month’s *Sargassum* biomass contributes to the nutrient pool mainly with a weight (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {W}_{{lag}} $$\end{document}Wlag) after part of it is lost minus a dilution contribution. We then used MATLAB’s fminsearch program to optimize the statistical model. To quantify uncertainty, the optimization was repeated independently 50 times. The mean weight distribution with error bars indicating the standard deviation across the 50 runs is shown in Supplementary Fig. 7. Result indicated that it would be possible to take the mean concentration over the previous 12 or 15 months as well, but using the three-months running mean concentrations around the 12-month lag simplifies the formulation while maintaining high skill metrics, making it easier to implement in a forecasting framework. #### Regional sensitivity We furthermore tested whenever the pre-conditioning of the GASB is linked mostly to *Sargassum* growing in the Tropical and Equatorial Atlantic region (Supplementary Fig. 8A). We used the ODATIS daily gridded *Sargassum* area coverage (10.12770/8fe1cdcb-f4ea-4c81-8543-50f0b39b4eca, accessed on 15-04-2026) derived from satellite imagery and extracted concentrations over the Tropical and Equatorial Atlantic (longitude from −55 to −20, latitude from 0 to 10) to re-run our statistical model modifying the self-fertilization term. The resulting timeseries compared to the observed ones with *R* = 0.76 and MSE = 3.57, indicating that this region explains a substantial portion of the *Sargassum* bloom observed across the broader area. This outcome can be expected because the area at any time accounts for about 50% of the total bloom (see for example newer datasets providing concentrations by subdivisions at [https://optics.marine.usf.edu/projects/saws.html](https://optics.marine.usf.edu/projects/saws.html)) and represents the area most impacted by upwelling (Fig. 2). Nonetheless R decreases and the fit deteriorates, especially in the last few years, when the impact of the Sargassumsphere is higher and the growth is considerable also in the region excluded from the calculations. We also tested using only fall and early winter conditions as bloom precursors, as suggested in ref. ^30^. In this case, the model performance degraded significantly (*R* = 0.55 and MSE = 5.85) (Supplementary Fig. 8B). Limiting the model to data from the GASB area alone - therefore limiting the region where MLD changes are calculated – caused only a non-significant deterioration (*R* = 0.84 and MSE = 4.37), maintaining a good fit especially from 2020 onward, supporting the hypothesis that the role of the Sargassumsphere has increased over time. ### Sensitivity to noise in the data To test the sensitivity to uncertainties in the datasets used in our regression model, we considered two of the most used reanalysis products and quantified their differences in the representation of SST and MLD with respect to GLORYS. The two ocean reanalysis products, ORAS5^40^ and SODA Version 3.3.1^41^, have comparable resolution being both created at 1/4° horizontal resolution, while GLORYS has a higher resolution of 1/12°. We considered eleven years, from January 1994 to December 2014, because this the common period to all three products in their consolidated (validated) version, and calculated the standard deviation of the deseasonalized and detrended time series of SST and MLD over the region considered in our work, finding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\sigma }_{{SST}}=0.04 $$\end{document}σSST=0.04 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\sigma }_{{MLD}}=1.2567 $$\end{document}σMLD=1.2567. We then performed MonteCarlo simulations accounting for this uncertainty (Supplementary Fig. 9). In relation to the *Sargassum* concentration data, we compared the ones we used, adapted from ref. ^30^, with those published in ref. ^1^, which however are limited to December 2018 (Supplementary Fig. 10). Their difference is small and depends on the retrieval algorithms and assumptions made therein, which are not provided in the referenced papers. We verified that limiting our model to 2018 and retraining it using each dataset, caused only a small change in the coefficients, without impacting the model skill. ### Properties of the regression model’s residuals We used the outcome of the MonteCarlo simulations to evaluate the residuals of the regression model, finding that they follow a Gaussian distribution as shown by the standard Quantile-Quantile (Q-Q) plot (Supplementary Fig. 11A). When considering single realizations, however, the variance of the residuals may show a dependence on the values of the independent variables, with extreme events (both very high and very low) producing the larger residuals and suggesting some degree of heteroscedasticity. In addition, the residual autocorrelation is large or moderate up to 2 months, indicating a weak positive linear dependence between the model’s error in any given period and its error 1-2 months earlier. Over longer lags or longer (we tested up to 16 months) the Ljung–Box Q test^45^ confirms that there is no statistically significant autocorrelation (*p* > 0.05). We stress that short-lag dependence in residuals is common in Earth system applications^46^, and some degree of heteroscedasticity and autocorrelation can reflect underlying structural changes and nonlinear responses in complex systems that may undergo transitions^47,48^, as in our case, and is not linked to a failure of the predictive model. ### Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. ## Supplementary information Supplementary Information Reporting Summary Transparent Peer Review file ## Source data Source Data