Optimal wave reflection as a mechanism for seagrass self-organization

Ecosystems threatened by climate change can boost their resilience by developing spatial patterns. Spatially regular patterns in wave-exposed seagrass meadows are attributed to self-organization, yet underlying mechanisms are not well understood. Here, we show that these patterns could emerge from feedbacks between wave reflection and seagrass-induced bedform growth. We derive a theoretical model for surface waves propagating over a growing seagrass bed. Wave-induced bed shear stress shapes bedforms which, in turn, trigger wave reflection. Numerical simulations show seagrass pattern development once wave forcing exceeds a critical amplitude. In line with Mediterranean Sea field observations, these patterns have half the wavelength of the forcing waves. Our results raise the hypothesis that pattern formation optimizes the potential of seagrass meadows to reflect wave energy, and a clear direction for future field campaigns. If wave-reflecting pattern formation increases ecosystem resilience under globally intensifying wave climates, these ecosystems may inspire nature-based coastal protection measures.

equation, and the dynamic boundary condition at the water surface is given by the Bernoulli equation, with g the gravitational acceleration.Here, atmospheric pressure was assumed to be constant and set equal to zero at the water surface, for convenience.The kinematic boundary condition at the water surface is given by with horizontal gradient operator ⃗ ∇ h = (∂ x , ∂ y ).Finally, the kinematic boundary condition at the bottom is given by We assume here that the location of the bottom boundary, z = −H, is linearly related to dimensionless seagrass density n(x, y,t), i.e.
where h is the constant water depth in the absence of seagrass (i.e., the location of the flat seabed) and s can be interpreted as a topography coefficient, translating seagrass density to topographic elevation.This parameter defines the rate or efficiency at which the seagrass meadow raises the seabed through the formation of rhizome-interwoven sediment deposits or "matte" 2 .Its value is the bed elevation gained due to the meadow when the (dimensionless) seagrass density is of order 1.Its exact value (see Supplementary Table S1) should be calibrated in follow-up studies, but is here chosen consistently with the observation that "mattes" can continue to develop up to more than 10m in height 3,4 .We furthermore assume that dimensionless seagrass density is governed by with net seagrass mortality rate ω (death rate minus growth rate), facilitative interaction coefficient α, competitive interaction coefficient β and lateral plant dispersion coefficient δ (all dimensionless).Equation (7) was introduced earlier 5 ; here we have ignored non-local interactions between plants and plant dispersion due to clonal growth.We assume that wave-induced bed shear stress τ b (x, y,t) exerted on the bottom boundary increases seagrass mortality, i.e.
where ω b is the constant background value of the net mortality rate and ω c represents the coupling strength between bed shear stress and seagrass mortality.Following earlier studies [6][7][8] , we express wave-induced bed shear stress as with (constant) water mass density ρ, horizontal near-bed orbital velocity amplitude b and wave friction factor f w .For simplicity, we assume f w is constant.Assuming that the horizontal orbital velocity is a sine wave with period T = 2π/σ (with angular wave frequency σ ), the horizontal velocity amplitude equals √ 2 times the root mean square of the horizontal orbital velocity.In that case, Then, the full equations for coupled water wave and seagrass dynamics are where D denotes the entire fluid domain, i.e. |x| < ∞, |y| < ∞, −H ≤ z ≤ η (Supplementary Figure S1), with H as in (6) and ω as in ( 8) - (10).Equations ( 11) -( 15) describe the spatio-temporal dynamics of the state vector ⃗ χ = {φ , η, n}.

Series expansion
Inspired by earlier studies 9 , we assume that system state ⃗ χ can be written as a series expansion, where ) and ε is a small and constant scaling parameter, |ε| << 1. Substitution of the above series expansion in full equations ( 11) - (15) yields O(1) equations for the basic state ⃗ χ 0 , i.e.
with O(1) seagrass mortality and bed shear stress Similarly, the O(ε) equations for perturbation state ⃗ χ 1 are with O(ε) bed shear stress and seagrass mortality In equation ( 26), the term s∂ t n 1 was neglected, given that the local rate of change in topographic elevation is typically much smaller than vertical water flow induced by topographic relief, ∂ z φ 1 .

Basic state solution
Under specific simplifying assumptions, the O(1) hydrodynamic equations ( 16) -(19) reduce to the linearized equations for surface gravity waves, which have a well-known analytical solution 1 .Firstly, the free surface at z = η can be replaced by a fixed surface at z = 0 under the assumptions that wave steepness is small and that the real water depth, η + H, is approximately equal to H, i.e.
with O(1) wave amplitude a and wavenumber κ.Secondly, requiring that equations ( 17) and ( 18) can be linearized and that ∂ t φ 0 is balanced by gη 0 yields two additional requirements, i.e.
We will first assume that the latter scaling relation is valid and will confirm this afterwards.Thirdly, we assume that n 0 is a constant and will verify this afterwards.Equations ( 16) -( 19) then can be linearized and written as where D 0 is as D but with vertical range −H 0 ≤ z ≤ 0. Here, bottom boundary level H could be approximated by a fixed reference level, since , by definition of the series expansion of n.
When assuming that φ 0 can be written as a plane wave with a fixed vertical structure, travelling in the horizontal direction, the well-known linear gravity wave solution is found, i.e.
This is a travelling plane harmonic wave with wave amplitude a, wavenumber κ and angular frequency σ .The latter equation is the dispersion relation.Note that the wave vector was chosen to point in the x-direction, without loss of generality.The first assumption in (31) can now be verified.The second requirement of (31) becomes which only poses an additional requirement with respect to (30) when κH 0 << 1.In this latter case (the shallow water limit), tanh(κH 0 ) ≈ κH 0 , such that a/H 0 << 1.In other words, the second requirement of (31) is consistent with (30).Under the assumption that the O(1) seagrass density is constant in space and time, equation (20) reduces to and the bed shear stress becomes (10) after O(1), we find Wavenumber κ and seagrass density n 0 can now be solved (numerically) from equations ( 39) and (41), after which state variables η 0 and φ 0 can also be found.The solution of the basic state, for one example of parameter settings, is shown in the Results section of the main manuscript.

Perturbation equations
Whereas the linearized equations derived in the previous subsection allow for determination of the linear stability of basic state ⃗ χ 0 , the linearisation assumptions inevitably break down when seagrass-induced topographic modulations continue to grow in amplitude.In the full nonlinear seagrass equation ( 15), the competitive interaction term −β n 3 ensures saturation of seagrass growth.We anticipate that, by allowing nonlinear terms in the linearized seagrass equation ( 27), saturation of seagrass-induced topographic modulations will ensure that the assumptions used to linearize the hydrodynamic equations ( 11) -( 14) remain valid.We expect that a set of fully linear perturbation-state wave equations combined with a nonlinear perturbation-state seagrass equation will lead to a physically appropriate quasi-equilibrium state, while avoiding the computational difficulties involved in numerically solving the fully nonlinear wave equations.
Thus, instead of assuming |n m | = O(ε m ) and expanding and we expect that |n 1 | ≲ |n 0 | due to the growth saturation term that will be introduced.Note that we have now directly assumed that n 0 is constant, as was found in the previous subsection.With this notation, the basic-and perturbation-state approximations of bottom boundary equation ( 14) become The term on the right-hand side of the basic-state equation can be neglected (consistent with equation 35), given that , which is much smaller than 1 since seagrass-induced topographic modulations have a typical amplitude |s n 1 | of about 1m, 10 .However, the term s ⃗ ∇ h φ 0 • ⃗ ∇ h n 1 does balance with ∂ z φ 1 in the perturbation equation, where it dominates over the term s In conclusion, we replace the fully linear perturbation equations (23 -27) by the following final perturbation equations, with linearized hydrodynamics and nonlinear seagrass dynamics: with τ b1 as in (28), ω 1 as in (29), and the perturbation-state facilitative interaction coefficient α 1 defined as Linear stability of the non-vegetated state Equation (41) has two physically relevant solutions, namely a non-vegetated state (n 0 = 0) and a vegetated state (n 0 > 0).As explained in the main manuscript, the linear stability analysis of the vegetated state is not trivial and is therefore inferred from numerical time-integration of perturbation equations (45 -49).These numerical analyses will be discussed in the next section.Linear stability of the non-vegetated state, however, can be calculated analytically, as will be shown here.Consider a small, uniform perturbation ν to any of the uniform equilibrium solutions n 0 of equation ( 41).The evolution of this uniform perturbation can be written as Using equations ( 41) and (42), we can write with When we assume that perturbing n 0 with ν does not lead to significant changes in κ, equation (51) can be written as 3 (54)

5/12
Expanding this expression around ν = 0 and truncating after first order yields the linear perturbation equation, i.e.
We define the critical wave amplitude a * as the value of a where the eigenvalue λ ν becomes zero, i.e.
For the current choice of real-valued parameters and ω b < 0 (see Supplementary Table S1), the critical wave amplitude is a * + and λ ν is positive for a < a * + and negative for a > a * + .That is, the non-vegetated solution (n 0 = 0) is linearly unstable when wave forcing amplitude is weaker than a * + and linearly stable when wave forcing is stronger than a * + .

Numerical solutions
In this section, we describe the methods used to numerically solve the uniform basic state ⃗ χ 0 and to numerically time-integrate the equations for perturbation state ⃗ χ 1 .

Numerical solution of the uniform basic state
To solve the uniform basic state, dispersion relation (39) and seagrass equation (41) need to be solved numerically for water wavenumber κ and seagrass density n 0 .Surface elevation η 0 and velocity potential φ 0 then follow from equations ( 37) and ( 38).
In this study we consider monochromatic waves, hence angular wave frequency σ is a constant.To fix σ , we solve dispersion relation (39) analytically by first assuming an absence of seagrass (i.e., H 0 = h) and choosing κn /2 as estimate for κ, where κn is the average seagrass pattern wavenumber observed around Mallorca 11 .Here, we have a priori hypothesized that the seagrass bedforms are formed by Bragg resonance (the hypothesis which is then tested in this study), such that the incoming water wavenumber is half the seagrass pattern wavenumber 9 .
We then numerically solve the uniform basic state as a function of wave amplitude a. Starting at a = 0, seagrass equation (41) can be solved analytically for n 0 .This value is then used to solve dispersion equation (39) numerically, using κn /2 as initial guess for κ.Then, wave amplitude a is increased in small steps.For each increment, the seagrass equation and dispersion relation are solved together, using the solution of n 0 and κ at the previous value of a as initial guess for the numerical solver.This procedure is used both for the unvegetated solution (n 0 = 0) and the vegetated solution (n 0 > 0).Parameters used to solve ⃗ χ 0 are listed in Supplementary Table S1, where [n] denotes the unit of seagrass density n (which, in our study, we assume to be dimensionless).
We obtained the values of the seagrass-related parameters (δ , ω b , ω c , α and β ) as follows.First, we started from the seagrass equation described in previous studies 5,11 , and simplified this equation by retaining only the terms responsible for facilitative and competitive interaction, local mortality and lateral diffusive dispersion, since we hypothesize that these terms alone (together with the effect of wave-induced bed shear stress on seagrass mortality) are sufficient to explain seagrass patterning.Second, we slightly adjusted the values of α and β , such that seagrass density has only one solution for each value of control parameter ω.This is done to avoid bistability, as was present in the original model 5 but which would distract the attention of our study from its focus on pattern formation.Third, we made the mortality rate ω explicitly dependent on the wave-induced bed shear stress, via equation (8).Herein, we chose a negative value for the background mortality, ω b , such that a vegetated meadow exists in the absence of wave forcing.The coupling strength ω c between bed shear stress and seagrass mortality was then chosen such that, for a realistic range of water wave amplitudes, the seagrass density gradually decreases and eventually becomes zero.Fourth, lateral seagrass dispersion coefficient δ was chosen such that the diffusion length scale, 2 δ /∆t, is of the order of the horizontal grid cell spacing ∆x, to avoid numerical instabilities.Finally, we assigned the units of seconds and meters to the independent variables time and space, respectively, in the seagrass equation, which was dimensionless in its original form 5 .Given this same time-and spatial scale, the seagrass dynamics can be directly coupled to the hydrodynamics.Although the current study is theoretical and aims purely at testing the hypothesis that wave reflection can explain seagrass self-organization, it is desirable to further calibrate the seagrass-related model parameters based on field or literature study in follow-up research.
Note that due to the modifications described above, the values of seagrass-related parameters (δ , ω b , ω c , etc.) are not directly comparable to the values reported in previous studies 5,11 .This is because, in our model, we need to bridge the gap between the hydrodynamic timescales (wave dynamics, on the order of seconds) and the biogeomorphic timescales (seagrass dynamics, on the order of months to years).We bridge this timescale difference by effectively speeding up the simulated (bio)geomorphic evolution.This approach is applied often in morphological modelling and consists in multiplying the (bio)geomorphic dynamics i.e. η 1 (x,t + ∆t) = η 1 (x,t) + ∆t • F η1 (x,t) at z = 0.However, to minimize reflection of the perturbation wave (φ 1 , η 1 ) on the lateral domain boundaries, a so-called sponge layer needs to be implemented first [15][16][17] .A sponge layer is implemented on both lateral ends of the domain to mimic an infinitely long domain.Within this sponge layer, a damping term ( f SL and f SR for the sponge layers on the left-and right-hand side, respectively) is imposed whose strength increases gradually from the interior of the domain towards the boundary, i.e.
with damping terms where A S is the amplitude of the damping functions and ξ L is the horizontal coordinate that linearly increases from 0 to 1 within the left sponge layer (i.e., increases in the negative x-direction) and ξ R linearly increases from 0 to 1 in the right sponge layer (i.e., increases in the positive x-direction), see Supplementary Figure S2d.Outside of the sponge layers, f SL and f SR are zero.P S is the power of the damping function, i.e. it determines how quickly the damping term increases within the sponge layer.Furthermore, η 1L and η 1R are the values of η 1 that are imposed at the outermost left and right boundary.We choose both values to be zero, such the perturbation wave field is completely damped out at the end of the sponge layers and wave reflection is minimized.Here, both sponge layers are chosen to fit exactly N S water wavelengths λ .Secondly, similarly to the kinematic surface boundary condition, the dynamic surface boundary condition, equation ( 46), is rewritten as and the same damping functions f SL and f SR are applied to dampen φ 1 (x, z = 0,t) towards a constant value of φ 1L = 1 and φ 1R = 1 at the outermost left and right boundaries, i.e.
The constant values φ 1L and φ 1R are imposed as boundary conditions on the outermost left and right boundaries, for all values of z.
Thirdly, the Laplace equation ( 45) is solved numerically for φ 1 .An iterative Laplacian inversion technique (Gauss-Seidel method with successive over-relaxation) is used 14 , where r is the residual used to let the iterative process converge and the relaxation parameter p r determines the speed of convergence.Laplacian inversion is performed for each time step.Initially, the solution of φ 1 from the previous time step is used as the first guess in the iterative inversion method.Then, the boundary conditions at the left and right boundaries (φ 1L and φ 1R ) are imposed as well as the surface boundary condition computed from the dynamic surface boundary condition (as discussed above).Then, one Laplacian inversion iteration is performed, only for the interior domain (i.e., excluding the left, right, surface and bottom boundaries).After that, the bottom boundary condition, equation (48), is discretized to obtain the bottom boundary condition for φ 1 .This boundary condition is then imposed as well.After this entire procedure, the relative error at iteration step k + 1 is calculated, i.e.
where φ 1 k+1 is the solution of the Laplacian inversion method after iteration k + 1.This solution then becomes the initial guess of φ 1 for the next iteration, and so on.This iteration loop is repeated until the relative error becomes smaller than a critical value e crit , i.e. until the solution has converged.
Fourthly, once the Laplacian inversion process has converged, the basic-state wave field (φ 0 , η 0 ) is updated (i.e., computed for this new time step, t + ∆t).Using these values and the updated values for φ 1 and η 1 , we calculate the instantaneous perturbation bed shear stress during this time step, which is proportional to the integrand in equation (28).
Fifthly, this entire process (i.e., numerically solving the perturbation wave field and calculating the basic-state wave field for a new time step) is repeated N ∆t times, i.e. for one entire wave period T .Hence, we now know the instantaneous perturbation bed shear stress values for all time steps within this entire wave period.Following equation (28), the perturbation bed shear stress averaged over this wave-period, τ b1 , can be computed.

) Note that the fact that |η 1 |
<< |η 0 | and |n 1 | ≲ |n 0 | justified approximating the water surface level z = η by a fixed level z = 0 and replacing the bottom boundary level z = −H by a constant level z = −H 0 in these perturbation equations.