A mathematical model for inter-speciﬁc interactions in seagrasses

Seagrasses are vital organisms in coastal waters, and the drastic demise of their population in the last decades has worrying implications for marine ecosystems. Spatial models for seagrass meadows provide a mathematical framework to study their dynamical processes and emergent collective behavior. These models are crucial to predict the response of seagrasses to diﬀerent global warming scenarios, analyze the resilience of existing seagrass distributions, and optimize restoration strategies. In this article, we propose a model that includes interactions among diﬀerent species based on the clonal growth of seagrasses. We present a theoretical analysis of the model considering the speciﬁc case of the seagrass-seaweed interaction between Cymodocea nodosa and Caulerpa prolifera . Our simulations successfully reproduce ﬁeld observations of shoot densities in mixed meadows in the Ebro River Delta in the Mediterranean Sea. Besides, the proposed model allows us to investigate the possible underlying mechanisms that mediate the interaction among the two macrophytes.


INTRODUCTION
Seagrass meadows are one of the most valuable structural elements in marine coastal areas since they provide ecosystemic services in the form of nutrient supply, refuge, and nursery ground to many species of fishes and invertebrates [Costanza et al., 1997].In addition to support marine biodiversity, they create architectural structure as benthic producers, contribute to the water quality and sediment stabilization [Orth et al., 2006], protect coasts lines from strong waves [Fonseca andCahalan, 1992, Sánchez-González et al., 2011], and are responsible for globally significant carbon sequestration [Duarte et al., 2005].Moreover, seagrass beds support fisheries and provide livelihoods for millions of people in coastal communities [Watson et al., 1993].Seagrass ecosystems, like many others, are under global threat due to anthropogenic impact [Orth et al., 2006,Hughes et al., 2009].Eutrophication, coastal development, water pollution, increased mooring activity, competition with invasive species, and global warming are some of the many causes that are leading seagrasses to experience an accelerating decline in the last decades [Waycott et al., 2009].
The critical demise of the worldwide population of seagrasses requires immediate and strategic actions for marine coastal protection.Mathematical models provide a theoretical framework to study the most relevant mechanisms that govern the dynamics of these ecosystems, and they can predict the possible future scenarios subjected to different environmental conditions.They can also evaluate the resilience of the existing population to stress factors and identify tipping points, where a slight change in conditions may cause irreversible loss.Thus, mathematical models constitute an essential tool to assess the health of ecosystems and generate more informed decisions for the sustainable management of marine coastal areas.An agent-based spatial model for seagrasses and other organisms that exhibit clonal growth was proposed in [Sintes et al., 2005].This model collects detailed information of each shoot, and by iterating a set of empirically-based clonal growth rules, reproduces the non-linearities in the dynamics of a meadow [Sintes et al., 2006].The applications of the numerical simulations produced by the model are diverse, and range from assessing the CO 2 capture potential by seagrasses [Duarte et al., 2013] to estimate the age of living meadows [Arnaud-Haond et al., 2012].A macroscopic description of the agent-based model was recently proposed to predict the spontaneous formation of spatial patterns in seagrass meadows at seascape level [Ruiz-Reynés et al., 2017].In this model the change in the seagrass density obeys a set of coupled first order partial differential equations that include non-linear local and non-local interaction terms.By losing the detailed information on the shoot and ramet development, this macroscopic model works efficiently at larger spatial scales and it can be used to study the underlying mechanisms behind self-organization [Ruiz-Reynés andGomila, 2019, Ruiz-Reynés et al., 2020].
Current models in the literature have not considered interactions among different species.Several interacting species commonly constitute seagrass ecosystems and form spatial distributions that range from perfectly separated domains to mixed meadows.The study of interspecific interactions is key to identify the main processes that shape the distribution of the seascape.For instance, the introduction of invasive species that compete for resources with native ones can drastically affect the habitat conditions [Norton, 1976, de Villèle and Verlaque, 1995, Alós et al., 2016].Also, the tropicalization of temperate latitudes due to global warming aggravates the spread of exotic species that are more resilient to higher water temperatures [Rius et al., 2014, Vergés et al., 2014, Wesselmann et al., 2021].Since each species have different responses to the thermal stress [Savva et al., 2018, Collier et al., 2011], the addition of inter-specific interactions to the present numerical model is necessary to predict the dynamics of seagrasses under different global warming scenarios.
In this work, we introduce a generalization of the agent-based model [Sintes et al., 2005] that includes a cross-interaction term among clonal species either seagrasses or seaweeds.Such interaction is implemented in the local interaction term through coupling parameters that quantify the strength of the interaction between any pair of species.These parameters cannot be directly measured and need to be inferred indirectly by field observations.We will test the model and explore the role of the coupling for the specific case of the seagrass-seaweed interaction between Cymodocea nodosa and Caulerpa prolifera.These two species commonly form mixed meadows, and several studies showed that they negatively influence each other [Tuya et al., 2013, Pérez-Ruzafa et al., 2012].In this article, we propose a systematic way to fix the interaction parameters of the model and reproduce field observations of mixed meadows of C. nodosa and C. prolifera in the Alfacs Bay (Ebro River Delta).

Numerical Model
In this article we propose a numerical model to study the interactions among different species of seagrasses.Our model follows similar clonal growth rules as those formulated in [Sintes et al., 2005, Sintes et al., 2006] for non-interacting seagrasses that we briefly summarized here.In their model, the development of clonal networks was simulated using a set of ecologically relevant parameters that can be easily derived from empirical observations, such as: the rhizome elongation rate [v], that sets the horizontal spread of the clone; the branching rate [ν 0 ], that controls the capacity of the clone to form dense networks; the branching angle [φ], that determines the efficiency of the space occupation; the spacer length [δ], that measures the length of the piece of rhizome between consecutive shoots; and the shoot mortality rate [µ] [Bell andTomlinson, 1980, Marbà andDuarte, 1998].The simulation starts placing a seed (a shoot carrying an apical meristem) and assigning to it a unitary vector û, randomly oriented, setting the direction of growth of the rhizome.At each iteration, the following steps are repeated: 1.A rhizome, that originates in a randomly selected apex, is proposed to grow from its current position, r 0 , to r = r 0 + δ û.The proposal is accepted if no other shoot is present within an exclusion area of radius σ < δ centered at r.The value of σ is set to avoid multiple shoots occupying the same position and to preserve the shoot density reported in natural stands of the species.Then, the apex is relocated at r where a new shoot will develop.In this process, the direction of growth, û, does not change.
2. Time is increased as ∆t = δ/(vN a (t)), where N a (t) is the number of apices at time t.
3. A new branch, holding a growing apex, may develop at r with a probability p ν (t) = ν 0 ∆tN a (t).
The new branch will extend a long a new unitary vector û forming an angle φ with û along the right or left side of û, randomly chosen.Only one branch is possible at the position where the apex is located.
4. Within this time step, a number of shoots are removed from the meadow with a probability p µ (t) = µ∆t/N s (t), where N s (t) is the number of shoots at time t.It is assumed the shoot mortality to be an age-independent event [Duarte et al., 1994].
In this work, an alternative and more efficient type of exclusion procedure is implemented.A characteristic value of the shoot density, according to the empirical observations, is set for each of the competing species, ρ max .A square grid, representing the transects placed in the meadow, is superimposed on top of the continuum space.The grid spacing is set to 20 cm.When a rhizome is proposed to extend into a cell which density has reached its saturation value ρ max , it will not advance in that iteration.This method generates the same results as the exclusion area principle and reduces substantially the computing time.
The application of the above mentioned growth rules lead to an homogeneous spatial development of the clones that strongly depend on the balance between the branching and the shoot mortality rates [Sintes et al., 2005], but is not able to reproduce the reported self-organized marine vegetation patterns, such as, stripes or fairy circles [Bonacorsi et al., 2013,Borum et al., 2013,Frederiksen et al., 2004, Pasqualini et al., 1999, Ruiz-Reynés et al., 2017, van der Heide et al., 2010].This problem has been addressed assuming non-linear density-dependent growth rates that include facilitation and competitive interactions [Ruiz-Reynés et al., 2017].Although the pattern formation in single specific meadows is directly related to the presence of non-local interactions, in the case of interacting species, and for the sake of simplicity, we will assume a local density dependent branching rate of the form: where ν 0 is the intrinsic branching rate that depends on external factors such as temperature or irradiance [Duarte, 1989].The local density dependence includes two terms: a facilitative one that is assumed to be linear and results from the positive contribution of the neighboring plants that dissipates the wave energy and prevents the shoot removal; and a non-linear competitive term that determines the environmental carrying capacity.The reduction in the branching rate may be the result of different mechanisms, such as competition for natural resources or self-shading in dense meadows [Invers et al., 1997].ρ = ρ/ρ max is the normalized local density, and α is a coefficient that controls the strength of the interaction.Given the parabolic shape of the interaction, the growth of over-and under-populated areas is penalized, whereas regions around an optimal density ρ = ρ max /2 is favoured.The shoot mortality rate, µ, is kept fixed.
Equation (1) can be easily extended to consider the interaction among N species as follows.We define a normalized local density for the species i = 1, . . ., N as: where ρ max,i is the saturation density for the i−species.γ ij is the coupling coefficient, which controls the strength of the competitive interaction between species i and j.The higher their value is, the more one species is affected by the presence of the others, and vice-versa.The case γ ij = 0, reduces to the single-species case.This normalized density ρi is used in Equation ( 1) to evaluate the branching rate for the i−species, ν i (ρ i ).The form of the interaction term is chosen to be simple yet generic, and its biological interpretation will depend on the specific species we consider.Once the system is initialized with a random distribution of seeds of the competing species, the simulation proceeds as follows: 1. Since species have different characteristic growing times: τ i = δ i /v i , at each iteration, one of the species is selected with probability 2. The rhizome that originates in the n th -apex of the i-species, randomly selected, is proposed to extend over a distance δ i û(n) i .
3. The apex will be relocated to its new position and a new shoot will develop only if the normalized local density in the corresponding cell, given by the Equation ( 2), fullfils: ρi < 1.

Time is increased by ∆t
, where N T a (t) is the total number of apices from all species at time t 5.A new branch with a growing apex will develop according to the branching rate , with N a,i (t) the number apices of the i−species.
6.During this time step, a number of shoots of the i−species are removed with probability p µ,i (t) = µ i ∆t/N s,i (t), with N s,i (t) the number shoots of the i−species.
The coupling coefficients γ i,j cannot be fixed by direct measurements and must be inferred indirectly by comparing the outcomes of the model to field observations.In this work, we will consider the seagrass-seaweed interaction between C. nodosa and C. prolifera and it will be used as a test case for the proposed model.The set of parameters used as model inputs for each of the selected species is presented the Supplementary Materials A, according to averaged experimental observations.The saturation density for both species is selected to be ρ max = 1800 shoots/m 2 .In the Supplementary Materials C, we have studied the response of the interacting model proposed above for single species, taking into account local self-interaction term, (2) with γ ij = 0. We find that the system tends to bi-stable density states that strongly depend on initial conditions.In the rest of this work, we will consider mortalities µ i < ν 0,i .This condition ensures that the species are out of the bi-stability region.
For two species (N = 2), the number of coupling parameters reduces to γ 12 and γ 21 .We will determine the values of these coefficients such that our simulations reproduce experimental observations.In our simulations we have used a system size of L × L, with L = 20 m and periodic boundary conditions.The results have been averaged over 15 independent realizations.

Experimental design
In order to evaluate the coupling parameters, the model outcome will be compared to the observed shoot densities in mixed meadows of C. nodosa and C. prolifera.The field observation was conducted in the Alfacs Bay (Ebro River Delta), a shallow embayment (up to 6 m deep) with an estimated surface area of 50 km 2 , located in the Northwestern Mediterranean.The fieldwork was performed on the northern shore of the embayment, which is covered by an extensive seagrass meadow of Cymodocea nodosa.The northern part of the bay receives seasonal freshwater inputs as runoff from rice paddy fields, with high nutrient and organic matter concentrations, and with suspended materials as well that increases water turbidity and, therefore, light deprivation, as compared to the southern part of the bay [Pérez et al., 1994].To avoid seasonal variability, plant sampling took place within a short period of time in two consecutive years (June 2018 and June 2019) when the seasonal peak in seagrass biomass and shoot density occurs [Mascaró et al., 2014].Points were randomly selected for a depth range, and distributed along, approximately, 10 km of seagrass meadow (See Supplementary Materials B for a map of the study site).At each sampling point, we collected with a hand-held corer (15 cm internal diameter) all the plant and algae biomass to a depth in the sediment of about 30 cm [Pérez et al., 2001].Each sample was rinsed in situ thoroughly with seawater to eliminate the sediment, and organic material sealed in plastic bags, immediately frozen and carried to the laboratory where it was conserved at −25 • C until processing.Samples from 2018 were transported to the laboratory.Laboratory procedures involved separating plant and algae material of each sample into the different fractions (shoots, roots and rhizomes for Cymodocea nodosa, and fronds, rhizoids and roots for Caulerpa prolifera) and counting the shoots and fronds, respectively.In 2019, in order to optimize sampling effort, the shoot/frond density was counted directly on board, and no samples were transported to the laboratory.

Symmetric interaction between two species
We consider a scenario in which two species, C. nodosa and C. prolifera, coexist in the same region.For simplicity, we start assuming symmetric couplings between them, i.e., γ 12 = γ 21 = γ.Our main result in this section is the identification of two distinct regions in the space of parameters: γ > 1, and γ < 1.From the definition of the coupling coefficients in (2), we observe that γ > 1 implies that the interaction between species is higher than the self-interaction.In this case, we expect meadows to separate into different domains to minimize the competition.The opposite happens for γ < 1, where the competition minimizes when the species mix.In the following, we study the outcomes of our mathematical model in these different regimes.All the simulations have been initialized with a random distribution of seeds that ensures a rather homogeneous spatial distribution.The initial shoot density for both species is ρ 0, Cn = ρ 0, Cp = 25 m −2 , where the subindexes Cn and Cp stand for C. nodosa and C. prolifera, respectively.The coefficient α has been set to α Cn = α Cp = 4.
The value of the shoot mortality rate for C. nodosa is fixed to its average field observation value µ Cn = 0.92 yr −1 (Supplementary Materials A), whereas the one for C. prolifera has been varied in the region µ Cp < ν 0,Cp in order investigate the change in the population of shoots and their spatial distribution.Here, we summarize the most relevant results of our simulations: γ > 1: In this case, the competition between shoots of different species is strengthened, and it favors the spatial segregation of shoots in different domains.The results for γ = 2 are shown in Figure 1.As expected, soon after the homogeneous initial condition has been set (t < 5 yr), well-defined domains separated by clear fronts emerge.For a shoot mortality rate for the C. prolifera, µ Cp > 0.73 yr −1 , the area covered by C. prolifera shrinks, and the space left is occupied by C. nodosa that will become the dominant species.The situation reverts for µ Cp < 0.7 yr −1 in which C. prolifera colonizes the whole meadow.The coexistence of different species is not possible and, after a transient period characterized by the spatial segregation of species, whose duration depends on µ Cp , the stable steady-state solution corresponds to a homogeneous mono-specific meadow.
γ < 1: The competition between species weakens, which favors the formation of mixed meadows.The results for γ = 0.5 are illustrated in Figure 2. The steady-state behaviour is quickly achieved (t < 2 yr) as seen in Figure 2 γ = 1: This is the limit case between the two previous regions (Figure 3).In this case, shoots of different species equally compete.The analysis of the saturation densities as a function of µ Cp (Figure 3(b)) shows that the coexistence region is extremely narrow in comparison with the one found for γ = 0.5 (Figure 2(b)).As a consequence, mixed meadows, that can last hundreds of years (see Figure 3

Experimental analysis and model validation
In this section, we study the response of the model to non-symmetric coupling coefficients γ 12 = γ 21 .One of our main observations is that the results from the symmetric case generalize and the formation of mixed meadows is guaranteed when the coupling coefficients are in the range 0 < γ ij < 1.In Fig. 4, we plot the shoot densities of C. nodosa (ρ Cn ) vs. C. prolifera (ρ Cp ) for different combinations of interaction parameters in 0 < γ ij < 1.We find that their steady-states are mixed meadow solutions that follow linear relations of the type ρ Cn = Aρ Cp +B, as in the symmetric case (Fig. 2(c)).In Fig. 4(a), we kept the mortality of C. nodosa (µ Cn ) fixed, while varying the mortality of C. prolifera (µ Cp ).Interestingly, we observe that the coupling coefficient γ 12 governs the magnitude of the slope A, which follows the relation γ 12 ∼ −A.We could expect this relation from the linear stability analysis of the model at equilibrium (See Supplementary Materials D).For small values of γ 12 , the density ρ Cn is less affected by a change in the mortality µ Cp , resulting in a much smaller slope A. Changing the other coefficient γ 21 , translates the points along the same line since a decrease in γ 21 favors the growth of C. prolifera at the expense of C. nodosa, and vice-versa.
An analogous analysis can be done for Fig. 4(b), where we kept the mortality µ Cp fixed and varied µ Cn .Here, the coefficient γ 21 controls the slope of the linear regression.In this case, it follows the inverse relation γ 21 ∼ −A −1 , and the coefficient γ 12 does not alter the slope.The different sets of points fit the linear regressions of slope A, with an error not higher than 2% in any of the cases.
After the analysis of the model behavior, it is possible to determine the coupling coefficients γ ij in real meadows by comparing the results of our simulations with the experimental data.We measured the averaged shoot density in meadows of C. nodosa and C. prolifera in 106 sampling points in the Alfacs Bay (Ebro River Delta -Spain), shown with purple dots in Figure 5.The presence of mixed meadows of these two species indicates values of the couplings γ 12 , γ 21 < 1.Although data is highly scattered, it is possible to find linear relationships between the shoot densities of C. nodosa, vs. C. prolifera: ρ Cn = Aρ Cp + B. We use the least-square method to fit the data [Watson, 1967].This method minimizes the errors of one of the data sets and considers the other as a control variable with no error.If we choose to minimize the errors in the measurements of ρ Cn , the best fit to the data gives a slope of A = −0.47 ± 0.11 (solid blue line).This slope, within its error, is in a very good agreement with the model outcome for γ 12 ∼ 0.5 according to Fig. 4(a).The linear regression that minimizes the errors in ρ Cp has a slope of A = −3.3± 0.7 (solid red line), which coincides with the simulations when γ 21 ∼ 0.3 (Fig. 4(b)).Therefore, the numerical model nicely reproduces field data for γ 12 = 0.5 and γ 21 = 0.3.For these coupling coefficients, the simulations with fixed µ Cn and varying µ Cp generate the dotted blue line, and similarly for the dotted red line, by keeping µ Cp constant and mo µ Cn .The match between model and data is one of our main results, which establishes γ ij as the key parameter to allow the connection between the numerical model and the experimental data and provide information on the strength of the interaction between species.In all previous studies we assumed an initial condition consisting of a mixed meadow with the same amount of seeds of both species homogeneously distributed all over the available space.We analyze the sensitivity of the results to the initial condition.We will assume the case in which patches of different species are initially placed apart from each other.In a square region of size L 2 = 40 × 40 m 2 , we locate a homogeneous density of seeds along two vertical stripes, 1 m wide, one of C. nodosa, centered at x = −10 m, and another of C. prolifera, centered at x = 10 m (Figure 6 (t = 0 yr)).During a short period of time, two separated mono-specific meadows develop.Shortly after that, the two domains collide forming a clear front (Figure 6 (t ∼ 8 yr)).At this stage, the two species have reached their maximum shoot density values.This front quickly dilutes giving rise to a mixed meadow of C. prolifera and C. nodosa (Figure 6 (t ∼ 19 yr)).The final state of the system is comparable to that collected assuming a homogeneous initial distribution of seeds (see Figure 2(d)).

DISCUSSION
In this paper we present the results of an agent-based model that investigates the dynamics of growth of clonal organisms (like seagrasses and seaweeds) in a scenario of species competition.The interaction among species has been implemented through a coupling parameter, γ ij , placed in the local interaction term (eq.( 2)) that quantifies the strength of the interaction between any pair of species.We found the value of this coupling parameter to determine the model outcome that ranges from separated and well defined mono-specific domains to mixed ones.A comparison between the model outcome and field measurements, conducted in the Alfacs Bay where large meadows of C. nodosa and C. prolifera are present, allowed us to determine the coupling coefficients between these two species (γ 12 = 0.5, γ 21 = 0.3) and it is one of the main findings of this paper (Fig. 5).In our simulations, the choice of coupling parameters γ ij is decisive since they reflect the type and magnitude of the interaction between species: γ 12 controls the influence of C. prolifera over C. nodosa, and vice-versa for γ 21 .The asymmetry found in the results is a clear indication that the negative effect of the presence of C. prolifera over C. nodosa is higher than vice-versa.The fact that both coefficients are in the γ ij < 1 regime implies that the seagrass/seaweed competes less with the shoots of the other species than with their own.Therefore, in order to minimize the competition, species self-organize and form mixed meadows.This effect should typically be explained by several biologically relevant mechanisms related to the anatomy, metabolism, or other properties of the plants.Although C. prolifera and C. nodosa have comparable magnitudes in their clonal growth rates (Supplementary Materials A), they have very different root lengths: while C. nodosa extends into the sediment to a mean depth of 14.1cm, the roots of C. prolifera just elongate an average of 3cm [Duarte et al., 1998, Bedinger et al., 2013].This characteristic allows the shoots of the two species to absorb the nutrients from different soil regions, making self-competition higher than the interspecific competition, favoring the appearance of mixed meadows.Therefore, our numerical analysis leads us to conjecture that the difference in root lengths between species is one of the main mechanisms that mediate the formation of mixed meadows of C. nodosa and C. prolifera.
Our analysis is consistent with previous results found in the literature.In [Tuya et al., 2013], shoot densities in mixed meadows of C. nodosa and C. prolifera were collected in 16 sites across Gran Canaria island (Spain, Atlantic Ocean) in the summer season, when the meadows should be at their optimal growth.The best linear fit to their data (ρ 1 = Aρ 2 + B), minimizing the errors in the measurements of C. nodosa, led to a slope A = −0.71± 0.35 which is consistent with our findings (Fig. 5).Despite the similarities found, different slopes can be expected due to the very different conditions between the locations of both experiments: the sites in our experiment were located in the Ebro River Delta, that is characterized by very calm, shallow, and high-on-nutrient waters, while sites in [Tuya et al., 2013] were in Gran Canary island, whose coastline has a lack of pronounced geographical accidents and it is very exposed to winds and currents.
In [Tuya et al., 2013], they also performed an experiment to establish the direct effect of C. prolifera over C. nodosa.In mixed meadows, 100% of the fronds of C. prolifera were removed.Eight months later, the density of shoots of C. nodosa increased 1.5-times.This result is also consistent with our findings.When the two species are allowed to grow separately, setting appropiate initial conditions as those depicted in Fig. 6, C. nodosa patches reach a maximum density of ρ Cn = 2000 m −2 (maximum of the solid blue line at t ∼ 8 yr).Once the two species meet and interact a stable mixed meadow of C. prolifera and C. nodosa develop with densities ρ Cn = 1480 m −2 and ρ Cp = 1270 m −2 .Thus, isolated meadows of C. nodosa have an average shoot density 1.35-times higher that in mixed meadows.Therefore, our model supports the observations made in [Tuya et al., 2013], where they conclude that the appearance of C. prolifera partially contributes to the demise of the population of C. nodosa.Also in [Pérez-Ruzafa et al., 2012], the changes of seagrass beds in Mar Menor Coastal Lagoon (Spain, Mediterranean Sea) were analyzed.They observed a decrease and almost total loss of C. nodosa in deeper areas of the lagoon (2 − 7m) after a colonization event of C. prolifera in the early 1970's.Our simulations sustain the hypothesis that the invasion of C. prolifera should negatively influence the abundance of C. nodosa.Still, this cannot solely explain the loss of C. nodosa, and we should consider other effects that deteriorate the seagrass.For example, light availability could be a key effect in the dynamics of C. nodosa and C. prolifera meadows, since it affects both species very differently.While C. prolifera suffers from photoinhibition in transparent shallow areas, the shoot density of C. nodosa could be limited by an increased content of silt, clay, organic matter, or phytoplankton [Pérez-Ruzafa et al., 2012, Pérez-Ruzafa et al., 2019, García-Sánchez et al., 2012].In [Belando et al., 2021], the most recent study of the distribution of C. nodosa and C. prolifera in the Mar Menor, they study the presence of mixed meadows in the deeper areas of the lagoon (less than −4m).Their data does not show a negative correlation between the biomass of both macrophytes, which is in disagreement with our findings in Fig. 5.The discrepancy between both observations might be related to the different environmental conditions among study sites, and it should be further investigated.
In this work, we have presented the first numerical model that studies local interactions among different species of rhizomatous macrophytes.Although the present work has focused on mixed meadows of C. nodosa and C. prolifera, the proposed model can be applied to study the interactions between generic N -species.For example, by properly adjusting the clonal growth parameters and coupling coefficients, we could simulate the dynamics of tropical seagrass beds, with a large diversity of coexisting species [Short et al., 2007].Also, this numerical model can be used to predict the re-distribution of seagrasses over the ocean floor due to global warming, paying special attention to the competition between species with different thermal thresholds, such as P. oceanica and C. nodosa, or the role of invasive species [Savva et al., 2018, Rius et al., 2014].Finally, this work sets the necessary framework to study long-range interactions between seagrass species.These types of interactions are often mediated by chemical concentrations and are responsible for the appearance of optimal sea-scape patterns [Ruiz-Reynés et al., 2017].In future work, we will extend the current model to consider such effects.
(a).In this regime, the averaged shoot densities for the different species and different values of the shoot mortality rate of C. prolifera, µ Cp , are shown in Figure 2(b).Interestingly, the saturation densities for the C. nodosa vs. C. prolifera in these mixed meadows follow a linear relationship (see Fig. 2(c)).The best fit to the data gives a slope of −0.48±0.01.This result indicates that the occupation of C. prolifera fronds increases by a factor of two as the shoots of C. nodosa decreases.A representative snapshot of a mixed meadow after t = 10 yr of growth at µ Cp = 0.65 yr −1 is shown in Figure 2(d).
(a)  for µ Cp = 0.70 yr −1 ), are a transient state towards a homogeneous mono-specific meadow.

Figure 1 :Figure 2 :Figure 3 :
Figure 1: Change in the spatial organization of C. nodosa and C. prolifera in a mixed meadow.The coupling coefficients are γ 12 = γ 21 = 2.The shoot mortality rate for C. nodosa is µ Cn = 0, 92 yr −1 , whereas the one for C. prolifera changes from µ Cp = 0.77 to 0.70 yr −1 , from top to bottom.Different snapshots are taken at t = 50, 300 yr.Regions colored in blue (red) are dominated by the presence of C. nodosa (C.prolifera ), and the color green represents regions of coexistence of both species.The color bar shows the difference between ρ Cp − ρ Cn , and it has units of shoots/m 2 .In the right column, the change in the average population of shoots for both species is shown.

Figure 4 :
Figure 4: Sensitivity analysis of the model for the coupling parameters: γ 12 , γ 21 < 1.We plot the shoot density of C. nodosa vs. C. prolifera in stable mixed meadows.Each color represents a different combination of the coupling coefficients.Each data point corresponds to a different value the C. prolifera mortality rate, µ Cp , in (a), and the C. nodosa mortality rate, µ Cn , in (b).The different sets of points fit the linear regressions of slope A, with an error not higher than 2% in any of the cases.

Figure 5 :
Figure 5: Comparison between the model outcome and the experimental data.The purple dots correspond to experimental measurements of shoot densities of C. nodosa and C. prolifera in mixed meadows in the Alfacs Bay (Ebro River Delta, Spain).The solid lines correspond to two types of linear fits of the experimental data, minimizing both the errors in the measurements of density of C. nodosa (solid blue line), or minimizing the errors in C. prolifera (solid red line).The dotted lines are the results of our simulations assuming coupling parameters γ 12 = 0.5, γ 21 = 0.3.The dotted blue line is reproduced by keeping the mortality of C.nodosa fixed, while varying the mortality C. prolifera, and vice-versa for the dotted red line.For a better comparison with the field measurements, the maximum shoot density has been adjusted to the average values found in the Alfacs Bay location (ρ max,Cn = 2300 m −2 and ρ max,Cp = 1900 m −2 ).

Figure 6 :
Figure 6: (Left) Spatial distribution, at different time steps, of two competing species: C. nodosa (blue) and C. prolifera (red), starting from two well defined striped patches separated from each other.As time evolves, both species cover all available space forming a mixed meadow (green).The color bar shows the difference between ρ Cp − ρ Cn , and it has units of shoots/m 2 .(Right) Evolution of the averaged density of both species.The interaction coefficients are set to γ 12 = 0.5 and γ 21 = 0.3, and the mortality rates, µ Cp = 0.65 yr −1 , and µ Cn = 0.92 yr −1 .