A NEMO-based model of Sargassum distribution in the tropical Atlantic: description of the model and sensitivity analysis (NEMO-Sarg1.0)

The tropical Atlantic has been facing a massive proliferation of Sargassum since 2011, with severe environmental and socioeconomic impacts. The development of large-scale modeling of Sargassum transport and physiology is essential to clarify the link between Sargassum distribution and environmental conditions, and to lay the groundwork for a seasonal forecast at the scale of the tropical Atlantic basin. We developed a modeling framework based on the Nucleus for European Modelling of the Ocean (NEMO) ocean model, which integrates transport by currents and waves, and physiology of Sargassum with varying internal nutrients quota, and considers stranding at the coast. The model is initialized from basin-scale satellite observations, and performance was assessed over the year 2017. Model parameters are calibrated through the analysis of a large ensemble of simulations, and the sensitivity to forcing fields like riverine nutrient inputs, atmospheric deposition, and waves is discussed. Overall, results demonstrate the ability of the model to reproduce and forecast the seasonal cycle and large-scale distribution of Sargassum biomass.

1 Introduction

The massive development of holopelagic Sargassum spp. in the northern tropical Atlantic Ocean from 2011 to the present has caused annual stranding in millions of tons on the coasts of the Lesser Antilles, Central America, Brazil, and western Africa (e.g., Smetacek and Zingone, 2013; Wang and Hu, 2016; Langin, 2018; Wang et al., 2019). The proliferation affects the whole tropical northern Atlantic area, as illustrated by satellite observations for summer 2017 (Fig. 1, Berline et al., 2020).

Figure 1 Sargassum fractional coverage obtained from MODIS in July–August 2017 (Berline et al., 2020, brown color scale; value between 0.001 % and 0.02 %) and surface chlorophyll distribution in July–August (green color scale; im mg m⁻³) based on GlobColour MODIS monthly product from 2010 to 2018. Circulation schematic of the surface currents is superimposed: the North Equatorial Current (NEC), the northern and southern branches of the South Equatorial Current (nSEC and sSEC), the North Equatorial Countercurrent (NECC), the North Brazil Current (NBC), the Caribbean Current (CC), and the Loop Current (LC).

Modeling and forecasting the Sargassum proliferation and strandings are essential for designing effective integrated risk management strategies and is a strong and pressing demand from the civil society. This operational challenge concerns both event forecasts (i.e., on a 1-week scale) and long-term forecasts (one to several months). While many efforts have been made for short-term forecasts, initiatives for reliable long-term forecasting are very scarce and face several scientific challenges such as the large uncertainties in Sargassum detection and biomass quantification (Wang et al., 2018; Ody et al., 2019), a lack of knowledge on Sargassum physiology, and last but not least the absence of tools specifically designed to reproduce the large-scale distribution of these macroalgae.

Recent studies suggest that the increasing incidence of Sargassum blooms and their year-to-year variability is multifactorial: it may result from riverine and atmospheric fertilization of the upper ocean, western Africa upwelling variability, vertical exchanges at the mixed-layer base in the region of the Intertropical Convergence Zone (ITCZ), or anomalous transport due to climate variability (Oviatt et al., 2019; Wang et al., 2019; Johns et al., 2020). This highlights the complexity of the phenomenon and the need for a basin-scale and interdisciplinary approach.

In the recent years, modeling effort mainly focused on the transport properties of Sargassum rafts by offshore currents (Wang and Hu, 2017; Brooks et al., 2018; Maréchal et al., 2017; Putman et al., 2018, 2020; Wang et al., 2019; Berline et al., 2020; Beron-Vera and Miron, 2020), with significant advances on the role of inertia in the drift trajectories (Brooks et al., 2019; Beron-Vera and Miron, 2020) and the importance of considering windage to properly resolve the drift of the Sargassum mats (Putman et al., 2020; Berline et al., 2020). To our knowledge, Brooks et al. (2018) were the first to integrate Sargassum physiology along the trajectories and showed that considering growth and mortality improved the modeling of the large-scale distribution of Sargassum. A similar result was obtained in Wang et al. (2019), although they did not consider directly the physiology of the algae but local growth rate based on satellite observations. Indeed, few studies have investigated the biology and ecology of this holopelagic Sargassum species that proliferate in the Atlantic and their response to the variability of environmental parameters (Lapointe, 1995, 1986; Hanisak and Samuel, 1987; Carpenter and Cox, 1974; Hanson, 1977; Howard and Menzies, 1969).

In the present paper, we describe the numerical model we developed to represent the distribution of holopelagic Sargassum. This model relies on an Eulerian approach and integrates both transport and a simplified physiology model of the macroalgae. It is based on the Nucleus for European Modelling of the Ocean (NEMO) modeling system, which is widely used by the research community and European ocean forecasting centers (e.g., Mercator Ocean International, ECMWF), allowing efficient parallelization and interfacing with physical–biogeochemical models. In the following section, we review current knowledge on the ecology of Sargassum. The modeling system is described in Sect. 3. Section 4 shows the performance of the model at seasonal scale and discusses sensitivity of the modeled Sargassum distribution to the forcing fields. Discussion and a summary are given in the final Section.

2 Physiological and ecological features of holopelagic Sargassum

Pelagic Sargassum species (to date Sargassum natans and S. fluitans) are brown algae (Phaeophyceae) that live at the surface of the ocean, never attached to any substrate. Within these two taxonomic groups, three types of Sargassum that can be distinguished according to morphological features appear to fuel the recent Sargassum inundations in the Caribbean: S. fluitans III, S. natans I, and S. natans VIII (Schell et al., 2015). We still lack knowledge on the distribution of these species, but in recent years, S. fluitans III was predominant in 2017 (Ody et al., 2019) and formed beaching on the Yucatán coast, comprising on average > 60 % of total wet biomass (García-Sánchez et al., 2020), whereas Schell et al. (2015) reported a predominance of S. natans VIII in 2015.

One individual Sargassum fragment can vary in length from just 1 cm to more than 1 m. Under the action of Langmuir cells and ocean currents, Sargassum tends to group together to form large floating rafts on the water surface (e.g., Langmuir, 1938; Zhong et al., 2012). Individuals in these aggregations can be easily dispersed when the dynamical conditions favorable to aggregation cease (Ody et al., 2019). These assemblages spread out horizontally and can reach several tens of kilometers and a few meters' thickness.

Biological and physiological features are species dependent. We know relatively little about the physiology of these Sargassum. Considering biomass, their maximum growth rate is estimated to be around 0.1 d⁻¹ (Lapointe, 1986; Hanisak and Samuel, 1987; Lapointe et al., 2014). The Sargassum growth is sensitive to light and temperature. Carpenter and Cox (1974) suggest light saturation under normal October light conditions in the Sargasso Sea (35 W m⁻²), while Hanisak and Samuel (1987) found a higher saturation range of ∼ 43–65 W m⁻². The temperature dependence in Hanisak and Samuel (1987) for Sargassum natans suggests a broad optimal temperature range of 18–30 ^∘C and indicates no growth at 12 ^∘C. We lack information on the Sargassum fluitans response to the variability of the environmental parameters.

Lapointe (1986) highlights a growth mainly limited by phosphate availability, while the presence of nitrifying epiphytes (Carpenter, 1972; Michotey et al., 2020) could be a non-negligible source of nitrogen for Sargassum, as could urea and ammonium excreted by fish (Lapointe et al., 2014). It is also likely that Sargassum are able to store some nutrients in their tissues, as do other brown algae (e.g., Hanisak, 1983). This hypothesis is supported by Lapointe (1995), whose measurements revealed variable elemental compositions between individuals sampled in neritic vs. oceanic waters. In addition, no macroherbivores control holopelagic development offshore by grazing (Butler et al., 1983).

3 The Sargassum modeling framework

Our modeling strategy relies on a physical–biogeochemical model that resolves currents and nutrient variability in the Atlantic. We choose to develop a regional configuration so the model can be tuned to the region specificity and can be used to perform sensitivity tests as discussed in Sect. 4. In addition, our approach is based on a Sargassum model that integrates transport, stranding, and physiology of the macroalgae in the ocean surface layer, forced with surface fields obtained from the physical–biogeochemical model.

The two models share the same horizontal domain and grid, and both are based on the NEMO modeling system version 4.0 (Madec and the NEMO team, 2016). They are not coupled assuming then that Sargassum does not compete with phytoplankton and heterotrophic bacteria for nutrient resources, and that they are not grazed by the herbivore compartments of the biogeochemical model.

3.1 The physical–biogeochemical model TATL025BIO

For the physical component of the simulation, we use the regional NEMO-based configuration described in Hernandez et al. (2016, 2017) and Radenac et al. (2020) that covers the tropical Atlantic between 35^∘ S and 35^∘ N and from 100^∘ W to 15^∘ E. The resolution of the horizontal grid is $\mathrm{1}/\mathrm{4}$^∘ and there are 75 vertical levels, 24 of which are in the upper 100 m of the ocean. The depth interval ranges from 1 m at the surface to about 10 m at 100 m depth. Interannual atmospheric fluxes of momentum, heat, and freshwater are derived from the DFS5.2 product (Dussin et al., 2016) using bulk formulae from Large and Yeager (2009). Temperature, salinity, currents, and sea level from the Mercator global reanalysis GLORYS2V4 (Storto et al., 2018) are used to force the model at the lateral boundaries. This configuration has proven to properly represent many aspects of the tropical Atlantic dynamics such as the Amazon plume extent (Hernandez et al., 2016), the large-scale circulation (Kounta et al., 2018) or the surface salinity variability (Awo et al., 2018).

The physical model is coupled to the PISCES (Pelagic Interaction Scheme for Carbon and Ecosystem Studies) biogeochemical model (Aumont et al., 2015) that simulates the biological production and the biogeochemical cycles of carbon, nitrogen, phosphorus, silica, and iron. We use the PISCES-Q version with variable stoichiometry described in Kwiatkowski et al. (2018), with an explicit representation of three phytoplankton size classes (picophytoplankton, nanophytoplankton, and microphytoplankton) and two zooplankton compartments (nanozooplankton and mesozooplankton). The model also includes three non-living compartments (dissolved organic matter and small and large sinking particles). The biogeochemical model is initialized and forced at the lateral boundaries with dissolved inorganic carbon, dissolved organic carbon, alkalinity, and iron obtained from stabilized climatological 3-D fields of the global standard configuration ORCA2 (Aumont and Bopp, 2006), and nitrate, phosphate, silicate, and dissolved oxygen from the World Ocean Atlas (WOA; Garcia et al., 2010) observation database. The model is run from 2006 to 2017, and daily physical and biogeochemical fields are extracted to force the Sargassum model.

Particular care has been given to the prescription of the atmospheric and riverine fluxes of nutrients. The river runoffs are based on daily fluxes from the ISBA-CTRIP reanalysis (Decharme et al., 2019), which has proven to accurately reproduce the interannual variability of the large rivers of the basin (e.g., see Giffard et al., 2019, for the Amazon River). The riverine nutrient fluxes concentrations are from the GLOBAL-NEWS2 dataset, corrected with in situ observations from the Amazon basin water resources observatory database (HYBAM, https://hybam.obs-mip.fr/, last access: 1 February 2020) for the Amazon, Orinoco, and Congo rivers. As in Aumont et al. (2015), we consider an atmospheric supply of P, Fe, and Si from dust deposition. Here, these fluxes are forced using monthly dry-plus-wet deposition products (DUDPWTSUM) from the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) data available on the NASA Giovanni website (http://disc.sci.gsfc.nasa.gov/giovanni, last access: 3 April 2020). Comparison with in situ observations of dust fluxes in Guyana and in Barbados leads to an excellent match with MERRA2 fluxes (Prospero et al., 2020). A climatological deposit of N is obtained from global climate simulations carried out with the ARPEGE-Climate model (Michou et al., 2020). An interactive aerosol scheme including nitrate and ammonium particles (Drugé et al., 2019) is included in ARPEGE-Climate, allowing us to produce fields of wet and dry deposition of nitrogen, ammonium, and ammonia. Figure 2 illustrates that dust and nitrogen fluxes to the ocean are strong in our region of interest, and most particularly in the ITCZ region where atmospheric convergence may focus the wet fluxes.

Figure 2 Dust fluxes to the ocean from MERRA2 reanalysis for the year 2017 (a) and nitrogen flux to the ocean from ARPEGE simulations averaged over the period 2010–2014. Nitrogen fluxes consider the dry and wet fluxes at the ocean surface of NO₃, NH₄, and NH₃. Dust fluxes consider the total dust dry-plus-wet deposition product (DUDPWTSUM) from MERRA2.

The modeled chlorophyll for the year 2017 is compared with GlobColour satellite estimates of chlorophyll for the same year (Fig. 3a, b). Model NO₃ and PO₄ concentrations for 2017 are compared with historical in situ measurements (Fig. 3c–f) from the GLODAPV2 database (Olsen et al., 2016). The model reproduces the major chlorophyll structures and in particular the contrast between the oligotrophic subtropical gyre and productive coastal and equatorial upwellings (Fig. 3a, b). As many other models (e.g., see the CMIP6 model evaluation by Séférian et al., 2020), it struggles to reproduce the offshore extent of the large river plumes and Guinea Dome productivity. Moreover, coastal upwellings tend to be too productive offshore or downstream (for the equatorial upwelling). But above all, it represents realistically the chlorophyll distribution in the region of the ITCZ (∼ 0–10^∘ N) and in the Caribbean Sea. As observed, nitrate concentrations are high in the upwelling areas (Fig. 3c, d) but weaker than observed off these regions. It is worth noticing that historical observations of surface nitrate concentrations in the tropical band show very heterogeneous and contrasted values between cruises, so the reliability of a nitrate climatology in this area remains uncertain (Fig. 3c). The model reproduces realistically the observed interhemispheric gradient of surface phosphate concentrations even though it likely overestimates areas of high phosphate content (Fig. 3e, f).

Figure 3 Annual mean surface chlorophyll concentrations (mg chl a m⁻³) from the GlobColour satellite product (a) and model (b) for the year 2017. Spatial distribution of surface NO₃ and PO₄ concentration (c, e; mmol m⁻³) from historical cruises of the GLODAPv2_2016 database (https://www.nodc.noaa.gov/ocads/oceans/GLODAPv2/, last access: 1 June 2020) and annual mean surface NO₃ and PO₄ distribution from the model (d, f).

3.2 The Sargassum model NEMO-Sarg1.0

The Sargassum model relies on the strategy used to represent the distribution of other macroalgae species and their transport by 2-D advection–diffusion equations (e.g., Martins and Marques, 2002; Solidoro et al., 1997; Perrot et al., 2014; Ren et al., 2014; Ménesguen et al., 2006; Bergamasco and Zago, 1999). Here, we also consider sink due to stranding at the coast. Growth is modeled as a function of internal reserves of nutrients, dissolved inorganic nutrients in the external medium, irradiance, and sea temperature. As the eco-phycological features are species dependent, the actual knowledge on the three morphotypes of holopelagic Sargassum does allow us to discriminate between them. Following the formalism given in Ren et al. (2014), the physiological behavior is described from three state variables: the content in carbon (C), nitrogen (N), and phosphorus (P), with local variations reflecting the difference between uptake and loss rates.

$$\begin{array}{c}{\displaystyle \frac{\partial C}{\partial t}}=U_{\mathrm{C}}-{\mathit{\varphi }}_{\mathrm{C}}\\ {\displaystyle \frac{\partial N}{\partial t}}=U_{\mathrm{N}}-{\mathit{\varphi }}_{\mathrm{N}}\\ {\displaystyle \frac{\partial P}{\partial t}}=U_{\mathrm{P}}-{\mathit{\varphi }}_{\mathrm{P}},\end{array}$$

where U_C, U_N and U_P are the uptake rates of carbon, nitrogen, and phosphorus, respectively, and Φ_C, Φ_N, Φ_P the loss rates.

The rate of carbon uptake reads as follows: $U_{\mathrm{C}}=C\cdot {\mathit{\mu }}_{max}\cdot f\left[T\right]\cdot f\left[I\right]\cdot f\left[Q_{\mathrm{N}}\right]\cdot f\left[Q_{\mathrm{P}}\right]$, with μ_max the maximum net carbon growth rate, and the four subsequent terms standing for uptake limitation by temperature (T), solar radiation (I), N quota (Q_N), and P quota (Q_P), respectively. N and P quotas represent the ratios of nitrogen and phosphorus to carbon in the organism and are computed as $\mathrm{N}/\mathrm{C}$ and $\mathrm{P}/\mathrm{C}$, respectively. The C content (C) can be directly converted to dry biomass considering a mean carbon-to-dry-weight ratio of 27 % (Wang et al., 2018).

The temperature dependence is adapted from Martins and Marques (2002):

$$f\left(T\right)=e^{-\frac{\mathrm{1}}{\mathrm{2}}\cdot {\left(\frac{T-T_{\mathrm{opt}}}{T_x-T}\right)}^{\mathrm{2}}},$$

with $T_x=T_{min}$ for T≤T_opt; $T_x=T_{max}$ for T>T_opt. T_opt is the optimum temperature at which growth rate is maximum, T_min is the lower temperature limit below which growth ceases, and T_max is the upper temperature limit above which growth ceases. Such a function aims at representing a broad optimal temperature range as suggested by experiments in Hanisak and Samuel (1987). The dependence on light is expressed as follows in order to mimic results from Hanisak and Samuel (1987):

$$f\left(I\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{\left(-\mathrm{0.1}\frac{I}{I_{\mathrm{opt}}}\right)}}}.$$

We have very little information on the response curve relating the nutrient quota to Sargassum growth but experiments for brown seaweed suggest a hyperbolic relationship (e.g., Hanisak, 1983). So, the dependence on the internal nitrogen and phosphorus pools is computed as a hyperbolic curve controlled by the minimum and maximum cell quotas:

$$\begin{array}{c}{\displaystyle }f\left(Q_{\mathrm{N}}\right)=\left({\displaystyle \frac{\mathrm{1}-Q_{\mathrm{Nmin}}/Q_{\mathrm{N}}}{\mathrm{1}-Q_{\mathrm{Nmin}}/Q_{\mathrm{Nmax}}}}\right)\\ {\displaystyle }f\left(Q_{\mathrm{P}}\right)=\left({\displaystyle \frac{\mathrm{1}-Q_{\mathrm{Pmin}}/Q_{\mathrm{P}}}{\mathrm{1}-Q_{\mathrm{Pmin}}/Q_{\mathrm{Pmax}}}}\right).\end{array}$$

The nitrogen and phosphorus uptake rates depend on the nitrogen (V_Nmax) and phosphorus (V_Pmax) maximum uptake velocities, a Monod kinetic that relates uptake to nutrient concentrations in the water, and a function of quota which aims at representing downregulation of the transport system for N and P when approaching the maximum quotas (Lehman et al., 1975):

$$\begin{array}{c}{\displaystyle }{U}_{\mathrm{N}}=V_{\mathrm{Nmax}}\cdot C\cdot \left({\displaystyle \frac{\left[\mathrm{N}\right]}{K_{\mathrm{N}}+\left[\mathrm{N}\right]}}\right)\cdot \left({\displaystyle \frac{Q_{\mathrm{Nmax}}-Q_{\mathrm{N}}}{Q_{\mathrm{Nmax}}-Q_{\mathrm{Nmin}}}}\right)\\ {\displaystyle }{U}_{\mathrm{P}}=V_{\mathrm{Pmax}}\cdot C\cdot \left({\displaystyle \frac{\left[\mathrm{P}\right]}{K_{\mathrm{P}}+\left[\mathrm{P}\right]}}\right)\cdot \left({\displaystyle \frac{Q_{\mathrm{Pmax}}-Q_{\mathrm{P}}}{Q_{\mathrm{Pmax}}-Q_{\mathrm{Pmin}}}}\right).\end{array}$$

The carbon loss aims at representing mortality, stranding, and sinking:

$$\begin{array}{rl}{\mathit{\varphi }}_{\mathrm{C}}& =C\cdot ({\displaystyle \frac{m}{e^{-{\mathit{\lambda }}_{\mathrm{m}}\cdot (T-\mathrm{30}{}^{\circ }\mathrm{C})}}}+{\mathit{\delta }}_{\mathrm{land}}\cdot {\mathit{\alpha }}_{\mathrm{grnd}}\\ & +{\displaystyle \frac{m_{\mathrm{LC}}}{e^{-{\mathit{\lambda }}_{\mathrm{mLC}}\cdot \left(H_{\mathrm{LC}}-\mathrm{100}m\right)}}}).\end{array}$$

The mortality term depends on m the mortality rate, λ_m a mortality coefficient, and temperature T. We thus represent thallus senescence and bacterial activity as a growing function of temperature (Bendoricchio et al., 1994; Ren et al., 2014).

The stranding is a function of α_grnd which is a rate of Sargassum stranding per unit of time, and δ_land which is defined as follows:

$$\left\{\begin{array}{ll}{\mathit{\delta }}_{\mathrm{land}}=\mathrm{1}& \text{if model grid cell is adjacent to two or}\\ & \text{more pixels of land},\\ {\mathit{\delta }}_{\mathrm{land}}=\mathrm{0}& \text{otherwise}.\end{array}\right.$$

The sinking rate of Sargassum is estimated as a function of the Langmuir cell length scale H_LC and a sinking coefficient m_LC. The aim is to reproduce possible Sargassum loss by Langmuir cell, as hypothesized by Johnson and Richardson (1977) or Woodcock (1993) to explain large amounts of Sargassum observed at the sea floor (Schoener and Rowe, 1970; Baker et al., 2018). Following, Axell et al. (2014), we estimate H_LC as the depth that a water parcel with kinetic energy $u_{\mathrm{s}}^{\mathrm{2}}/\mathrm{2}$ can reach on its own by converting its kinetic energy to potential energy. This corresponds to $-{\int }_{-H_{\mathrm{LC}}}^{\mathrm{0}}{N}^{\mathrm{2}}\left(z\right)\cdot z\cdot \mathrm{d}z=\frac{\mathrm{1}}{\mathrm{2}}\cdot u_{\mathrm{s}}^{\mathrm{2}}$, with u_s the Stokes drift and N² the Brunt–Vaisala frequency. We fixed a 100 m depth threshold as Sargassum becomes massively negatively buoyant at these depths (Johnson and Richardson, 1977).

Losses of nitrate and phosphate are function of the loss of biomass and internal N and P quotas:

$${\mathit{\varphi }}_{\mathrm{N}}={\mathit{\varphi }}_{\mathrm{C}}\cdot Q_{\mathrm{N}},{\mathit{\varphi }}_{\mathrm{P}}={\mathit{\varphi }}_{\mathrm{C}}\cdot Q_{\mathrm{P}}.$$

The transport of C, N, and P is resolved using 2-D advection–diffusion equations discretized on a grid at $\mathrm{1}/\mathrm{4}$^∘ resolution with a single vertical layer representing a surface layer of water of 1 m depth. The surface velocities used for the transport account for surface currents, windage effect, and wave transport by Stokes drift:

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}{\mathit{\varphi }}_{\mathrm{transport}}\left(\mathrm{Nutrient}\right)& =-U\cdot {\displaystyle \frac{\partial \mathrm{Nutrient}}{\partial x}}-V\cdot {\displaystyle \frac{\partial \mathrm{Nutrient}}{\partial y}}\\ & +K_h\cdot {\mathrm{\nabla }}_h^{\mathrm{2}}\mathrm{Nutrient},\mathrm{with}\end{array}\\ {\displaystyle }\begin{array}{rl}(U,V)& =(u_{\mathrm{NEMO}},v_{\mathrm{NEMO}})+{\mathit{\alpha }}_{\mathrm{win}}\cdot (u_{\mathrm{10}\mathrm{m}},v_{\mathrm{10}\mathrm{m}})\\ & +(u_{\mathrm{Stokes}},v_{\mathrm{Stokes}}),\end{array}\end{array}$$

where (u_NEMO,v_NEMO) are the horizontal velocity obtained from the physical–biogeochemical model, α_win is a windage coefficient, (u_10 m,v_10 m) the components of the wind field at 10 m above the sea level, (u_Stokes,v_Stokes) the Stokes velocity, and K_h a diffusion coefficient.

3.3 Optimization and sensitivity experiments

The model simulations are performed for the year 2017 because basin-scale Sargassum fractional coverage observations from MODIS were available (Berline et al., 2020), with concurrent observations carried out during two cruises in the tropical Atlantic (Ody et al., 2019).

3.3.1 Initialization and forcing

The simulations are initialized using January Sargassum mean fractional coverage, converted into dry weight biomass considering a surface density of 3.34 kg m⁻² and then into carbon content C considering a mean carbon-to-dry-weight ratio of 27 % (Wang et al., 2018). The initial N and P content in Sargassum is derived from the initial C content and N and P quotas computed as the averaged values between their respective minimum values (Q_Nmin, Q_Nmin) and maximum values (Q_Nmax, Q_Nmax). The transport is forced by daily velocities from TATL025BIO simulations (see Sect. 3.1). The windage is forced with 3 h winds from the DFS5.2 dataset (Dussin et al., 2016), which were used to force the physical model. The Stokes drift in the surface layer is computed according to Breivik et al. (2014) and forced with hourly ERA5 Stokes drift product. Daily temperature, available irradiation, and Langmuir depth were also obtained from TATL025BIO. The seawater concentrations in [N] and [P] were obtained from TATL025BIO as the sum of NO₃ and NH₄ for [N], and PO₄ for [P], in the top surface layer.

Table 1 Sargassum model parameters.

3.3.2 Ensemble strategy

The Sargassum model is controlled by a large number (n=18) of physiological and physical parameters for which large uncertainties exist or most often have not been measured for the Sargassum species considered here. An ensemble approach has been adopted to adjust the set of parameters. We produced 10 000 sets of parameters with uniform distribution obtained from Latin hypercube sampling with multi-dimensional uniformity (Deutsch and Deutsch, 2012). These sets of parameters are generated on ranges of values obtained from the literature, when available (Table 1).

Figure 4 Seasonal distribution of Sargassum fractional coverage for the year 2017 from observations (left) and from a selected ensemble of 100 simulations (Ω₁₀₀) with varying parameters (right).

The range of maximum growth rate is derived from Lapointe et al. (2014) who observed maximum growth rates of Sargassum fluitans and Sargassum natans in neritic waters between 0.03 and 0.09 doubling d⁻¹. The set of parameters for temperature limitation is mainly derived from the results of Hanisak and Samuel (1987). In particular, the T_max range (40–50 ^∘C) is chosen to have a slight decrease of the limitation term at T>T_opt as observed in Hanisak and Samuel (1987). The N and P quotas are based on the observations by Lapointe (1995) from which we can estimate that on average $\mathrm{C}/\mathrm{N}$ ratios vary between 20 and 30 in neritic waters and between 40 and 70 in oceanic waters, while $\mathrm{C}/\mathrm{P}$ ratios vary between 200 and 500 in neritic waters and 700 and 1000 in oceanic waters. The lower and upper limits for the maximum nitrate uptake rate ([$\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$] mg N (mg C)⁻¹ d⁻¹) are estimated from measurements by Lapointe (1995). From this study, we estimate maximum carbon uptake rates in neritic water at ∼ 2 mg C (g dry wt)⁻¹ h⁻¹ with $\mathrm{C}/\mathrm{N}$ ratio of 20 and maximum carbon uptake rates in oceanic water at ∼ 1 mg C (g dry wt)⁻¹ h⁻¹ with $\mathrm{C}/\mathrm{N}$ ratio of 60. Since the model does not take into account a diurnal cycle of light, the maximum uptake has been divided by 3 in order to have a daily-mean maximum uptake. This ratio of 3 is obtained by comparing instantaneous gross production at full irradiance vs. measured gross production from culturing the Sargassum under natural irradiance in Lapointe et al. (2014). Similarly, the lower and upper limits for the maximum phosphorus uptake rate ($[\mathrm{9}\times {\mathrm{10}}^{-\mathrm{5}},\mathrm{7}\times {\mathrm{10}}^{-\mathrm{4}}]$ mg N (mg C)⁻¹ d⁻¹) are estimated from measurements by Lapointe (1995), which gives a $\mathrm{C}/\mathrm{P}$ ratio of 200 in neritic waters and 800 in oceanic waters.

Figure 5 Seasonal distribution of Sargassum fractional coverage for the year 2017 from observations (blue) and from a selected ensemble of 100 simulations (Ω₁₀₀) with varying parameters (mean in red and variance in shaded red) averaged in different areas: (a) tropical Atlantic [0–30^∘ N, 98–10^∘ W], (b) Caribbean Sea [8–22^∘ N, 85–55^∘ W], (c) Sargasso Sea [23–30^∘ N, 80–50^∘ W], and (d) eastern tropical North Atlantic [0–15^∘ N, 30^∘ W–0^∘ E].

3.3.3 Likelihood function

The likelihood function (L) to be minimized is the centered root mean square error between monthly series of observed and modeled Sargassum biomass contained in the tropical North Atlantic Ocean, defined as [0–30^∘ N, 98–10^∘ W], hereinafter NTAO:

$$L={\left[{\displaystyle \frac{\mathrm{1}}{N_t}}{\sum }_{t=\mathrm{1}}^{N_t}{\left[\left({\mathit{\psi }}_{\mathrm{m}}\left(t\right)-\stackrel{\mathrm{‾}}{{\mathit{\psi }}_{\mathrm{m}}}\right)-({\mathit{\psi }}_{\mathrm{obs}}\left(t\right)-\stackrel{\mathrm{‾}}{{\mathit{\psi }}_{\mathrm{obs}}})\right]}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}},$$

with N_t the number of outputs (12 months in our case), ψ_m(t) and ψ_obs(t) the modeled and observed biomass in the NTAO, and the overbars designating the annual average.

The choice of this likelihood function exerts no constraints on the spatial distribution of the Sargassum and does not consider uncertainties in the satellite measurements due to false detection, cloud masking, or Sargassum immersion. But as shown in the following section, such a simple strategy allows us to efficiently select a set of parameters that allow a good representation of the seasonal Sargassum distribution.

3.3.4 Sensitivity analysis

Once a set of values minimizing L was found, one-at-a-time sensitivity experiments were also performed, where only one single parameter is varied by ±10 % in two model runs while the others are fixed (e.g., Ren et al., 2014; Perrot et al., 2014). This allows us to capture the direct contribution from each parameter to the output variance, with parameters varying within an acceptable range. The set of values for the fixed parameters is given in Table 1 and were taken from the simulation with higher L. Following Ren et al. (2014), deviation from the baseline simulation (D) is quantified as

$$D={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{t=\mathrm{1}}^{t=\mathrm{12}}{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{m}(i,t)}-{\mathit{\psi }}_{\mathrm{m}\left(t\right)}^{\mathrm{0}}}{{\mathit{\psi }}_{\mathrm{m}\left(t\right)}^{\mathrm{0}}}},\text{with}N=\mathrm{12},$$

where ${\mathit{\psi }}_{\mathrm{m}\left(t\right)}^{\mathrm{0}}$ refers to the simulated biomass in the NTAO from the baseline parameter set at month t, and ψ_m(i,t) is the simulated biomass with one perturbed parameter i at month t. Two model runs were conducted with ±10 % change to the baseline value.

Figure 6 Likelihood (L) as a function of the different parameters that have been varied in the optimization experiment, shown here for the ensemble Ω₁₀₀ that composes the ensemble averages in Figs. 4 and 5. Units are given in Table 1.

4 Results

4.1 Seasonal distribution and sensitivity to model parameters

The observed and model seasonal distributions of Sargassum in 2017 are shown in Fig. 4. Model distribution is obtained from the ensemble mean of the 100 simulations with the lowest L. Hereafter, this ensemble will be referred to as Ω₁₀₀. Averages over selected areas are shown in Fig. 5. At this stage, it is worth recording that the selection of the ensemble simulations is performed without constraints on the spatial distribution, the only constraint being on the basin-scale seasonal biomass average (L).

Initialized in January, the model reproduces the seasonal distributions of Sargassum fairly closely. It simulates the Sargassum drift toward the Caribbean Sea, with a summer peak of biomass followed by a decrease until the end of the year, although the increase of biomass comes a little too early in February–March. In the Caribbean Sea, the largest abundance of Sargassum in the northern part of the basin near the Greater Antilles compared to the south of the Caribbean is consistent with observations. Despite a bloom that is occurring too early (May–June) near the Lesser Antilles, the modeled seasonal cycle is also consistent with the observations in this area (Fig. 5b). This is encouraging from the perspective of predicting strandings on the Caribbean islands. It also succeeds in maintaining the biomass in the Eastern part of the tropical North Atlantic below the ITCZ (Figs. 4 and 5d). In the Sargasso Sea area (Fig. 5c), simulations and observations consistently show an increase at the end of the year. The model tends to reproduce heavy proliferations in March–June which seem not to be observed. Given current knowledge, it is difficult to determine the causes of such a bias. It could be due to a bias in the nutrient content simulated by PISCES-Q at this period. Moreover, error in the Sargassum initial conditions (January) and in the transport parameterization can lead to this production too far north during March–June. An observation bias cannot be ruled out either since this area is very cloudy and presents very contrasted Sargassum aggregation properties.

The distribution of parameters for the ensemble Ω₁₀₀ is given in Fig. 6. There is a significant dispersion of most of the parameters, suggesting either low sensitivity to the parameters in question (as discussed below) or interdependency between them. The analysis shows that there are very different sets of parameters, within the prescribed ranges, that lead to similar seasonal biomass distribution (Fig. 5). Having said that, we observe that the mortality parameters m and λ_m minimizing L are biased toward low and high values, respectively. This highlights the key importance of this mortality function in representing the seasonal distribution. The half-saturation constants (K_N and K_P) and the maximum uptake rates of nutrients (V_Nmax and V_Pmax) are also distributed toward low and high values, respectively, suggesting some sensitivity of the Sargassum distribution to nutrient resources. These parameters are poorly constrained by observations, and such exercise will allow us to refine the optimization ranges for future studies.

The relative mean deviation of the Sargassum seasonal biomass to 10 % variations in model parameters shown in Fig. 7 confirms findings from the analysis of parameter dispersion in Fig. 6. The most influential parameters are the growth rate, the optimal irradiance, mortality dependence on temperature, and parameter controlling the nitrogen uptake V_Nmax. It also highlights the influence of the minimum N quota (Q_Nmin). There is moderate dependence on temperature which is in good agreement with measurements from Hanisak and Samuel (1987). In agreement with previous Lagrangian studies (Berline et al., 2020; Putman et al., 2020), we find some sensitivity to the windage parameter.

Figure 7 Sensitivity analysis to model parameters expressed as the mean relative deviation (%) between the baseline simulation and the simulation in which one parameter was modified by ±10 %. The set of parameters for the baseline simulation is given in Table 1.

The simulated stranding from the ensemble Ω₁₀₀ is shown in Fig. 8. At the scale of the tropical Atlantic, these strandings sum to 1.05 million tons of dry weight for the whole year. The predicted strandings in the Caribbean, northern Brazil, French Guiana, and Sierra Leone correspond to our current knowledge of Sargassum invasions (Bernard et al., 2019; Louime et al., 2017; Oviatt et al., 2019; Sissini et al., 2017; Oyesiku and Egunyomi, 2014; Smetacek and Zingone, 2013). Nevertheless, there are no available large-scale coastal observations or estimates of stranding to go further in the validation of these simulated strandings.

Figure 8 Cumulated annual Sargassum wet biomass stranding per area of 25×25 km for the year 2017 obtained from the ensemble Ω₁₀₀.

A peculiarity of our modeling strategy is to consider the Stokes drift. The Stokes drift induces a displacement of material parallel to the direction of wave propagation which directly transports the Sargassum. An ensemble of 100 simulations without Stokes drift (“NoStokes”), considering the set of physiological parameters from Ω₁₀₀, has therefore been conducted. Anomalies of annual Sargassum distribution with respect to Ω₁₀₀ are shown in Fig. 9a. Sargassum coverage is significantly increased in the central Atlantic but decreases sharply in the Caribbean and the southwestern part of the domain. This highlights the influence of waves, most probably due to the trade winds, in shaping the seasonal distribution, and transporting the algae southward.

Figure 9 Ensemble anomalies of Sargassum coverage (in % of surface) from sensitivity ensemble simulations of 100 members each, which were performed with the set of parameters from Ω₁₀₀.

4.2 Sensitivity to external nutrient forcing

We now use the Sargassum model to explore how and to which extent continental nutrient sources (riverine nutrient fluxes, dust deposition, atmospheric N deposition) could participate in the proliferation and may shape the seasonal distribution of Sargassum. First, TATL025BIO simulations were run by deactivating the river sources (simulation “noriver”), atmospheric deposition of P, Si, and Fe (simulation “nodust”), or atmospheric deposition of N (simulations “noNdepo”). We choose to distinguish between dust and N deposition because they do not have the same origin and because spatial, seasonal, and long-term variability is not necessarily the same. The simulations were produced from 2006, so the long-term influence of these forcings on the biogeochemical content is considered there. These simulations were then used to force three ensembles of 100 simulations which used the sets of parameters from Ω₁₀₀.

As expected for the noriver ensemble, the tropical western Atlantic, which is under influence of the Amazon and Orinoco rivers, is experiencing a decrease in nutrient concentrations (Fig. 10b). Nutrients also show a decrease in the equatorial Atlantic. But for other regions, such as the Sargasso Sea or the Guinea Dome, the long-term equilibration results in an increase in nitrogen concentration. The resulting Sargassum coverage shows a negative anomaly in the Caribbean Sea, in the Gulf of Mexico, and in the region of the ITCZ, with an annual mean basin-scale biomass decrease (Fig. 9a). For this ensemble, strandings are decreased from 0.65 to 0.56 Mt (−13 %).

Figure 10 Biogeochemical response to sensitivity experiments to river, dust, and N deposition: surface anomalies of NO₃+NH₄ (a–c) and PO₄ (d–f) with respect to the reference TATL025BIO simulation.

The nodust ensemble Sargassum distribution shows a slight positive anomaly over most of the domain (Fig. 9c), particularly in the central Atlantic and off western Africa, with a slight decrease in the Caribbean Sea and the Gulf of Mexico. Cumulative strandings sum to 0.72 Mt, and so, are slightly increased compared to the baseline ensemble (+10 %). This sensitivity can be explained by the large-scale phosphorus increase and regional increases in nitrogen and phosphorus concentrations (Fig. 10) in the biogeochemical simulation. This sensitivity is expected to be produced by the reduced iron concentrations, which limit the phytoplankton growth and thus the nutrient uptake.

The removal of atmospheric nitrogen deposition (noNdepo) leads to a global decrease in surface nitrogen concentration and an increase in surface phosphorus concentrations across the entire domain (Fig. 10). The resulting Sargassum coverage is significantly decreased over the whole domain, and the annual averaged stranding decreases by 1 % (Fig. 9d).

5 Discussion and summary

Since 2011, unprecedented massive strandings of the holopelagic Sargassum have been reported on the coasts of the Caribbean Sea, northern Brazil, and western Africa. In this paper, we developed an Eulerian model of Sargassum, which integrates transport, strandings and algal physiology. The Sargassum model is based on the ocean modeling platform NEMO and is forced by the physical and biogeochemical fields of a regional model (TATL025BIO), as well as by the ERA5 wave and wind fields. An ensemble approach has been used to optimize the physiological parameters. The results demonstrate the ability of the model to represent the spatial distribution and seasonal cycle of the Sargassum biomass in the western Atlantic and the Caribbean Sea.

While windage and inertial effects are considered of importance for the drift properties and large-scale advection (Brooks et al., 2019; Berline et al., 2020; Beron-Vera et al., 2020; Putman et al., 2018, 2020), we show here that Stokes drift has also significant impacts on the distribution of the Sargassum and in particular on their entrance in the Caribbean Sea. In addition to the anomalous currents that may be at the origin of the Sargassum bloom in 2011 (Johns et al., 2020), wave drift could also have contributed to the dissemination of the algae toward the equatorial Atlantic in the early 2010s. Wave transport of algae is therefore an important component of Sargassum modeling that has not yet been accounted for in previous modeling efforts (Brooks et al., 2018; Putman et al., 2018; Wang et al., 2019; Johns et al., 2020) and should deserve further attention.

Transport properties may also be impacted by the numerical choices and model resolution. Our model resolution is intermediate (∼ eddy permitting), so we lack some energy at the mesoscale. Since this mesoscale is particularly important for the dynamics in the Caribbean, Gulf of Mexico, or the North Brazil Current area, we would expect more realistic transport properties at higher resolution. But our experience is that $\mathrm{1}/\mathrm{4}$^∘ NEMO simulations work well in the region on many aspects of the regional dynamics, such as river plume extent (Hernandez et al., 2016, 2017), large-scale currents (Kounta et al., 2018), biogeochemistry (Radenac et al., 2020), and large-scale salinity distribution (Awo et al., 2018), among others. One reason is that the scales of variability in the tropics are larger than at midlatitudes. This is a posteriori confirmed by the present study since we show that the simulated ocean dynamics are good enough to represent the accumulation of Sargassum in the ITCZ, the advection in the Caribbean through the Antilles, and the episodic shedding of Loop Current eddies in the Gulf of Mexico. We also expect that model resolution is only part of the story regarding the dependence of the transport properties to numerics. Surface transport also depends on the vertical resolution of the model in the mixed layer, the vertical mixing scheme, the degree of coupling of the ocean circulation with the atmosphere or the waves, the wind product used to force the model, etc. In our model, the windage transport coefficient acts as an empirical factor that compensates lacking the explicit simulation of some of these processes and probably helps us to properly simulate a realistic large-scale Sargassum advection. Overall, we definitely need to rely on dedicated Lagrangian studies such as the one performed by Putman et al. (2018), Putman and He (2013), Berline et al. (2020), and Putman et al. (2020) to better constrain our model, and learn about best practices in terms of forcing Sargassum transport.

Figure 11 Average N and P uptake by phytoplankton (in mmol m⁻² d⁻¹; a and b, respectively) and by Sargassum (in µmol m⁻² d⁻¹; c and d, respectively). The N and P uptake by phytoplankton were obtained from the biogeochemical simulation, assuming a constant stoichiometry. The uptake rates were integrated over the model mixed-layer depth for each month of 2017 and averaged over the year. The bottom row shows the Sargassum vs. phytoplankton mean consumption ratio of (e) N and (f) P (in ‰) for the year 2017.

The ability of the model to simulate the large-scale distribution was also used to conduct sensitivity tests on the nutrient forcing from rivers and dust and atmospheric deposition. Here, it is worth remembering that the Sargassum model is not coupled with the biogeochemical model so it is not directly forced by these external inputs of nutrients but through the biogeochemical model. This prevents the representation of some opportunistic utilization of nutrients that could be done by the algae. Moreover, it is worth mentioning that the $\mathrm{N}/\mathrm{P}$ half-saturation constants obtained from the basin-scale optimization procedure are low (likely because the biogeochemical model tends to have low surface nutrient concentrations in the tropical northern Atlantic). This could limit the sensitivity of the model to high nutrient inputs. With these limitations in mind, we found a 17 % and 21 % decrease in annual Sargassum distribution in the experiments without river nutrient runoff and without atmospheric nitrogen deposition, respectively. This suggest that these forcings alone cannot fuel the total Sargassum biomass. Regarding the nutrient brought by the Amazon, this is in agreement with recent conclusions from Johns et al. (2020) and Jouanno et al. (2021), who suggest that the riverine fertilization of the tropical Atlantic is not at the origin of the phenomenon nor control its year-to-year variability. At this stage, the processes controlling the interannual variability and overall increase of Sargassum remains an open question that will deserve further attention. Application of the numerical tracer method initially proposed by Ménesguen et al. (2006), which tracks nitrogen or phosphorus from any source throughout the biogeochemical network, could help identify the nutrient sources that control the phenomenon without altering the large-scale biogeochemical content.

Several aspects which could be of potential importance for Sargassum growth have not been considered here. First, growth and mortality could depend on the age of the fragments, through colonization by epiphytes. There is a lack of knowledge on these aspects, and mesocosm experiments would be useful to better constrain such dependence in the model, if relevant. Second, we assume that Sargassum does not compete with phytoplankton for resources. The annual mean consumption of N and P in the mixed layer is shown in Fig. 11 for both phytoplankton and Sargassum, obtained from the NEMO-Sarg1.0 and PISCES-Q models, respectively. It reveals that the consumption of N and P by the phytoplankton is 2–3 orders of magnitude larger than the consumption of N and P by the Sargassum. So, at the basin scale, it seems reasonable to consider that Sargassum growth does not affect phytoplankton growth. But at the local scale (scale of a raft or scale of a bay) and particularly with low mixing conditions, Sargassum could compete with phytoplankton for resources. A full coupling of our Sargassum model with PISCES-Q may allow us to address such questions. Third, our results show a strong dependence on the nitrogen uptake parameters. We do not consider possible fixation of atmospheric N through diazotrophic assemblage. Biological nitrogen fixation by diazotrophic macroalgal association has been shown to be important for some Sargassum species (e.g. Sargassum horneri, Raut et al., 2018) and this could also be the case for the holopelagic Sargassum where epibionts N-fixating bacteria have been observed on these species (Carpenter, 1972; Michotey et al., 2020).

Finally, our modeling system succeeds in maintaining some biomass in the tropical central and eastern Atlantic. This pool of Sargassum has been shown to be of key importance in the year-to-year maintenance of the population (Wang et al., 2019). So, we expect the model to be useful to address the question of interannual variations.