Seagrass deformation affects fluid instability and tracer exchange in canopy flow

Monami is the synchronous waving of a submerged seagrass bed in response to unidirectional fluid flow. Here we develop a multiphase model for the dynamical instabilities and flow-driven collective motions of buoyant, deformable seagrass. We show that the impedance to flow due to the seagrass results in an unstable velocity shear layer at the canopy interface, leading to a periodic array of vortices that propagate downstream. Our simplified model, configured for unidirectional flow in a channel, provides a better understanding of the interaction between these vortices and the seagrass bed. Each passing vortex locally weakens the along-stream velocity at the canopy top, reducing the drag and allowing the deformed grass to straighten up just beneath it. This causes the grass to oscillate periodically even in the absence of water waves. Crucially, the maximal grass deflection is out of phase with the vortices. A phase diagram for the onset of instability shows its dependence on the fluid Reynolds number and an effective buoyancy parameter. Less buoyant grass is more easily deformed by the flow and forms a weaker shear layer, with smaller vortices and less material exchange across the canopy top. While higher Reynolds number leads to stronger vortices and larger waving amplitudes of the seagrass, waving amplitude is maximized at intermediate grass buoyancy. All together, our theory and computations develop an updated schematic of the instability mechanism consistent with experimental observations.


Introduction
Seagrass is typically deformable, which allows the grass blades to reconfigure according to the fluid load [Vogel, 2020].While emergent canopies -those that are in the inter-tidal zone and emerge above the water surface -need stiffness for the stems to stand up out of the water, fully submerged seagrass species (such as Halodule wrightii and Syringodium filiforme) tend to stand up by buoyancy [Wilson et al., 2010].In order to photosynthesize, submerged canopies have a typical height comparable to the water depth [Marion et al., 2014], which results in a significant portion of the flow being obstructed by the canopy.Seagrass beds exhibit a particularly rich set of dynamic behaviors due to their collective interaction with the flow.Hydrodynamic processes resulting from these interactions influence environmental processes such as sedimentation, transport of dissolved oxygen [Long et al., 2020] and nutrients, plant growth, and biomass production [Fonseca Figure 1: Schematic of the domain used in the simulations.The steady-state horizontal velocity profile u(z) is imposed as the velocity inlet boundary condition.A conformal map accounts for variations in h and h g to separate the overflow (where F = 0, in blue) and seagrass (where F = 0, in green) regions of the domain.Buoyant grass blades deform by the flow, apply a drag on the fluid F, and the composite tip positions determine h g .and Kenworthy, 1987;Grizzle et al., 1996;Nepf, 1999Nepf, , 2012]].Seagrass meadows are also believed to influence sediment deposition and resuspension [Short and Short, 1984;Walker et al., 1996], as vegetation can trap suspended materials [Short and Short, 1984] and reduce sediment movement [Fonseca and Fisher, 1986].
Instabilities of flow through submerged canopies yield a phenomenon known as monami -the progressive, synchronous oscillation of aquatic vegetation [Ackerman and Okubo, 1993;Nepf, 2012].Current explanations of monami [Ikeda and Kanazawa, 1996;Raupach et al., 1996;Ghisalberti and Nepf, 2002] rely on the existence of a shear layer at the top of the grass bed due to vegetation drag.Through a mechanism similar to the Kelvin-Helmholtz instability [Singh et al., 2016], the enhanced velocity shear near the grass top creates a sheet of vorticity that destabilizes into vortices over time.These vortices perturb the flow, which locally changes the deformation of grass blades and leads to synchronous oscillations of the grass bed.These perturbations to the mean flow have been observed experimentally and feature sweeps and ejections that occur at the leading and trailing edges of vortices, respectively [Ghisalberti and Nepf, 2006;Nezu and Sanjou, 2008;Okamoto and Nezu, 2009].Transport of material across the canopy has also been studied experimentally [Ghisalberti and Nepf, 2005;Nepf and Ghisalberti, 2008].Our numerical simulations provide a complementary and comprehensive picture of the fluid instability, vortex-seagrass interaction and tracer exchange between the seagrass bed and the overflow, in terms of its dependence on seagrass buoyancy and Reynolds number.
There are numerous modeling challenges in capturing the properties of this system, primarily related to the feedback mechanism between flow and vegetation.In a two-way coupled dynamic model, the fluid will apply a load on each vegetative structure, which causes a resultant deformation that, in turn, affects the flow [de Langre, 2008].Thus, in general, the fluid flow must be solved simultaneously with the configuration of each structure.These challenges have demanded sophisticated studies, both experimental [Dunn et al., 1996;Ghisalberti and Nepf, 2004;Okamoto and Nezu, 2009;Hu et al., 2014;Mandel et al., 2019] and numerical [Dupont et al., 2010;Zeller et al., 2015;Beudin et al., 2017;Mattis et al., 2019;Sundin and Bagheri, 2019].Most previous simplified models fall into one of two categories: models of flow over a specified set of rigid obstacles [Ghisalberti and Nepf, 2004;Singh et al., 2016], or models where grass deformation can occur, but does not change the flow profile [Luhar and Nepf, 2011].Fewer models [Wong et al., 2020] emphasize the coupling between grass deformation and flow, and our study is unique as it presents numerical simulations of the coupled system and uses them to study monami.
This study builds on previous work by Singh et al. [2016] and Wong et al. [2020] in analyzing the dynamics of flow through a submerged seagrass canopy and its resultant instabilities.Although monami is manifested in the grass motion, the drag exerted by the vegetation on the flow is central to the instability, and the resulting flow structures persist in laboratory experiments even when deformable grass is replaced by rigid dowels [Ghisalberti andNepf, 2002, 2006].Singh et al. [2016] proposed the seagrass effect on the fluid to be modeled as a continuum drag acting perpendicular to the blade, proportional to the number of stems per unit area, and established the dependence between viscous effects and flow instabilities by performing a linear stability analysis of flows through an array of rigid beams.Wong et al. [2020] expanded this model to account for flexible beams, derived the coupled equations of motion and relevant dimensionless groups, and performed a stability analysis to investigate conditions for the onset of instabilities.
We are interested in the impact that the grass blade deflection has on the onset of the instability, progression of the developed vortices, and material transport resulting from this interaction.Our model incorporates blade deformability into the two-phase model by Singh et al. [2016], but as opposed to the approach adopted by Wong et al. [2020], where the grass blades are modeled as linearly elastic flexible beams with one end clamped perpendicularly to the seabed, in our model the submerged grass blades stand up by buoyancy, do not resist shear (zero flexural rigidity), and are always in equilibrium with the flow (no contribution of the inertial term in the equations of motion for the grass).These assumptions simplify the equations of motion for the grass, while successfully reproducing the monami dynamics.
To model a submerged seagrass bed, we solve the Navier-Stokes equations for two phases: the grass-free overflow, and the grass-bed in which the seagrass contributes a bulk volumetric drag F that depends on the blade positions and velocity field (Fig. 1).The drag is quadratic in the velocity normal to the grass blades and hence depends on the grass shape, which in turn depends on the fluid drag.We model the shape of representative grass blades rooted to the bed in the center of each grid cell column (in plan view) by assuming a balance between drag, which deforms the grass blade, and buoyancy, which restores its shape to vertical.There are N grass blades per unit area that impose drag on the fluid, but do not block the flow.
Simulations are performed using a version of the non-hydrostatic Process Study Ocean Model (PSOM) [Mahadevan et al., 1996a,b].The submerged seagrass bed of undisturbed height is modeled in an open channel of undisturbed water height H using a grid that conforms to the free surface h(x, t) and seagrass height h g (x, t) as seen in Fig. 1.The along-channel coordinate is x, the vertical coordinate is z, and for the study described here, variations in cross-channel (y) direction are set to zero.The inflow velocity profile is in equilibrium with the grass, and within a buffer of the outflow boundary, we restore the velocity profile to the same equilibrium profile.All variables are non-dimensionalized using the undisturbed water height H as the characteristic length scale, the horizontal flow speed at the free surface U as the velocity scale, and H/U as the timescale.Variables are henceforth presented in dimensionless form.The dimensionless parameters that govern the solution are These are similar to Wong et al. [2020], except for β.Here, ν * is the constant eddy viscosity, ρ is the fluid density, ρ g is the grass density, g is acceleration due to gravity, c D is the quadratic drag coefficient, d is the thickness of the grass blades in the along-flow (x) direction, while b is the width of the grass blades (in the y-direction).The Reynolds (Re) and Froude (F r) numbers are standard parameters.The height ratio (r) is chosen as 0.5 in all our simulations, and this does not play a dominant role in any of the overall observations presented in this paper.The parameter λ governs the drag impedance by the seagrass and affects the velocity shear and is chosen as 1.The buoyancy parameter β is the ratio of seagrass buoyancy to drag and influences the shape of the grass blades.More buoyant grass has larger β and deforms less due to the flow.Choosing b and d as independent parameters enables us to change β without affecting λ.In this study, we perform numerical experiments for a range of Re, which is varied by changing ν * , and a range of β, which is varied by changing d.We analyze the onset of flow instability and amplitude of grass oscillations as a function of Reynolds number and grass deformation.We examine how the vortex position aligns with the shape of the grass bed and how the exchange of material across the grass canopy is affected by these dynamics.

Instability onset and progression
Before the instability onset, the flow and grass are both steady.The steady-state solution is a function of z alone, and can be calculated with a simplified one-dimensional coupled model that eliminates dependence in x and t (Methods section).The steady-state velocity profile u(z), grass shape (x g , z g ), and corresponding blade angle with the vertical θ(z), calculated for Re = 1000, r = 0.5, λ = 1, F r 2 = 0.1, and a range of values of the buoyancy parameter β, which in our model quantifies to what extent the blade can deform, is presented in Fig. 2. Solutions are computed with fluid boundary conditions u = 0 (no-slip) at the bottom (z = 0) and du/dz = 0 at the surface (z = 1).The horizontal pressure gradient is adjusted so that u = 1 at the surface.For the grass, the tension is zero at the tip, and the position is fixed at the bottom.The overbar is used to represent the steady-state solutions, which are independent of x and t.
Whether the grass relies on bending stiffness (as in Wong et al. [2020]) or buoyancy (this study) to restore its shape, the shape of the deformed grass and its implication for flow instability and fluid exchange are qualitatively similar (Supplementary information).The smaller β, the more blade deflection, the larger the angle θ along the blade, and the smaller the steady-state height h g corresponding to the height of the tip.As β increases, the velocity shear du/dz at the canopy top (z = h g ) monotonically grows, with the limiting case β → ∞ corresponding to a fixed, vertical blade, and maximum shear.
We initialize the channel model with the steady-state solution described above, and with no vertical velocity.For sufficiently large Re and λ, the shear layer at the canopy top is unstable, and instabilities are triggered spontaneously after finite time.The simulation (Fig. 3) with Re = 1000, r = 0.5, λ = 1, F r 2 = 0.1, and β = 0.10 exhibits the instability onset at t = 50 as seen in the vorticity ζ = ∂ z u − ∂ x w and vertical velocity w fields (Figs 3(a, c)).In most of the domain, w ≈ 0 and ζ ≈ du/dz with maximum ζ right above the canopy top, but some oscillations in ζ for x ∈ [16,30], where shear-instabilities start to grow and induce alternating vertical velocities.
When the instability is fully developed at t = 500 (Figs 3(b, d)), the vorticity rolls up to form vortices.Vortices are shed from x ≈ 7, grow until x ≈ 16, and stabilize in size as they propagate downstream with the flow (time evolution video in Supplementary information).The vortex centers lie between maxima and minima of vertical velocity, which peak just above the canopy (Fig. 3(d)).We also observe alternating vorticity maxima and minima near the seabed at z = 0 (Fig. 3(b)), which indicates that flow perturbations induced by the vortices penetrate to the bottom and cause flow reversal (u < 0) near the seabed, where the unperturbed velocity is already small due to the no-slip bottom boundary condition.
The space-time evolution of vertical velocity at the canopy top z = h g , where h g is the level corresponding to the steady-state canopy height (Fig. 4(a)), reveals an initial instability onset that originates around x ≈ 20 and t ≈ 50 that is swept out of the domain and replaced by unperturbed flow (video in Supplementary information).Another instability starts at t ≈ 100, now closer to the inlet at x ≈ 10, and develops to generate vortices almost periodically.Other than small variations that occur over time, notably a weakened vertical velocity field at x ≈ 30, t ≈ 300 and some fluctuation in the onset position, vortices are shed at a regular frequency proportional to the local fluid speed divided by the momentum shear layer thickness [Ho and Huerre, 1984] similar to Kelvin-Helmholtz instability.The vortices lead to alternating positive and negative w signatures that propagate downstream.
The generated vortices have constant speed of propagation (0.6), period (4.5), and wavelength (2.7).The propagation speed corresponds the mean downstream velocity at the height of the vortex centers.The The instability onset and strength of vortices can be assessed via the domain-wide rms of the vertical velocity , where N i and N k are the number of grid cells in x and z, respectively.The vertical velocity rms not only tracks when an instability has developed, but also indicates the strength of the vortices.The seagrass bed's response to the instability is assessed by the rms value of the grass height perturbation with respect to the steady state height h g over all blade representatives, δh rms g and it quantifies the vertical amplitude of grass blade oscillations.
The two quantities w rms and δh rms g , plotted as a function of t in Fig. 4(e), are zero before the onset, when there is no vertical velocity and all blades are at the steady-state shape.Their values increase rapidly at the initial onset of instability at t = 50, they decrease as the initial instability is swept from the domain, then increase again and assume nearly steady values for t > 150 as the long-term instability sets in.Because of the observed plateauing behavior of both curves, long-term values for w rms and δh rms g are defined as the time-average for t ∈ [300, 500] and used to inter-compare different cases in a parametric study.

Monami kinematics: effect of vortices on seagrass
We analyze the flow field, grass deflections, and free-surface height to evaluate how the instability interacts with the seagrass meadow to produce the oscillatory motion known as monami.The shear-driven instability induces a velocity perturbation field u = u − u with respect to the steady-state flow u = (u(z), 0) that deflects the grass blades from their steady-state position.
We subsample the domain to visualize two vortices at t = 500 in Fig. 5(a−c).We find that h(x) and h g (x) are approximately sinusoidal and out of phase, with peaks of h g slightly lagging troughs of h (Fig. 5(a)).While h has a more symmetric, sine-like profile, h g is less symmetric, with a steeper increase than decrease.The vortex cores are centered below the troughs of h and above the peaks of h g (Fig. 5(a, b)).The grass height h g is shown in Fig. 5(b) as a thick black line, and we observe that the clockwise vortices induce the grass blades beneath to straighten up.Near the seabed, directly below the vortices, negative vorticity values (in blue) indicate that horizontal velocity perturbations are strong enough to reverse the direction of the flow near the no-slip bottom.This generates convergence (and divergence) sites along the bed resulting in flow separation points that propagate with the vortices and could export sediment from the seabed.The velocity perturbation field calculated with respect to u(z) (Fig. 5(c)) further helps to visualize the response of grass blades to the flow perturbation.In Fig. 5(c), blades are colored based on the angle of deviation from the steady-state angle along the blade, δθ = θ − θ.There are two clockwise eddies that appear in the velocity perturbation field that align with the high vorticity regions in Fig. 5(b).Additionally, there are counter-clockwise vortices in between the vortex roll up highlighted by the velocity perturbation field, that induce a forward deflection of the blades.Blades immediately below the clockwise-vortices straighten up, while blades in between those vortices are subject to the action of counterclockwise-vortices that induce more deformation.
The distribution of h g (Fig. 5(d)) spans [0.475, 0.495] and is asymmetric with respect to the steady-state h g , showing that forward and downward deflection is more common and stronger than upward deflection.This is due to the fact that the drag force that deflects the grass acts normal to the blade.When the grass is downward deflected, the downward vertical velocity helps to enhance the deflection.When the grass is upward deflected with respect to its steady-state shape, the upward velocities are more or less parallel to the grass and do not contribute as much to the grass deflection as the downward velocity.The perturbation in the velocity field does not ever reverse the flow in the upper part of the canopy, and as a result the grass blades never move left from the vertical position.A schematic of how the vortices induce seagrass motion is presented in Fig. 5(e).The increased downward deflection of the grass occurs ahead and behind the vortex, where the counter-clockwise perturbation to the mean flow and the downward velocity cause a greater drag on the blades and deflect them forward and downward from their steady-state position (Fig. 5(c)).We identify the sweeps and ejections as corresponding to the perturbed velocity field immediately ahead and behind the vortices, around the seagrass height level (Fig. 5(b)).In the sweep region (ahead), the stronger velocity has a downward component and increases downward deflection of the grass.In the ejection region (behind), the weaker velocity has an upward component and induces an upward deflection of the grass.

Dependence of instability onset and waving amplitude on Reynolds number and grass buoyancy
The shear-driven instability and monami occur only when the drag-induced shear is strong enough for the vortex sheet at the canopy top to become unstable.This occurs when the velocity is large enough, i.e. above some critical value of Re, and when the grass-induced drag (λ) is sufficiently large.However, we find that β, the buoyancy parameter, also affects the instability onset and size of vortices.We use the long-term w rms as an indicator of instability and vortex strength, and δh rms g as an indicator of seagrass waving, to conduct a parametric study in which we vary Re and β over a range of values: Re ∈ [500, 1500] and β ∈ [0.02, 0.20].We run a total of 110 simulations, keeping the other parameters constant: r = 0.5, λ = 1, and F r 2 = 0.1.
We find no instability (w rms = 0) for small Re and β and that w rms increases for increasing values of Re and β, saturating for high values of both parameters (Fig. 6(a)).Larger w rms corresponds to stronger vortices inducing vertical velocities of greater magnitude, and it is reasonable to expect that as Re increases for fixed β the shear layer becomes stronger as do the resulting vortices.As β is increased (for fixed Re), the more buoyant blades are less deflected from the vertical.This sharpens the mean velocity gradient ( du/dz ) and results in stronger instabilities.Less buoyant blades are more deflected for the same Re, impose less drag (which is largely due to the velocity component normal to the grass) and inhibit the development of instabilities for fixed Re.Though our model uses buoyancy, the result is that the greater the deflection of stems, the less strong the instability, regardless of whether the deflection results from weak bending stiffness or buoyancy.
The critical combinations (Re, β) in the instability diagram above which w rms > 0 define an instability curve in Fig. 6(a).This curve is in agreement with the result that shear at the top of the canopy is the relevant criterion in determining the stability of steady unidirectional flows [Wong et al., 2020], as the velocity shear magnitude grows with both Re and β.
While w rms grows monotonically with Re and β (Fig. 6(a)), a non-monotonic behavior of the amplitude of waving, assessed by δh rms g , is observed in Fig. 6(b).In general, the amplitude of grass motion or δh rms g increases with Re.However, the maximum δh rms g occurs for intermediate values of β.Small β suppresses the shear instability and creates small w rms , whereas for larger values of β (e.g.0.20), the buoyancy of the grass resists its deformation, even though the fluid instability and w rms are stronger.For large β, the grass is almost vertical and the vortices do not induce an observable oscillatory motion.Grass oscillations are therefore maximized for specific combinations of (Re, β), with maxima δh rms g observed for high Re with intermediate β values.The observed trends point to the fact that experiments using rigid dowels [Ghisalberti andNepf, 2002, 2006] may overestimate the strength of the induced vortices, or not accurately predict their development, compared to what would be observed for more realistic, deformable seagrass beds.

Material exchange across the grass bed
To evaluate the impact of the vortices and the grass deflection on material exchange between the seagrass bed and the overflow, we model a tracer field and evaluate its transport.The initial tracer distribution C 0 (z) = 2 h g − z is linear in z with C 0 = 0 at the canopy top z = h g , C 0 > 0 within the grass bed and C 0 < 0 in the overflow.
The tracer transport and exchange is assessed for the same range of Re and β as above.A snapshot of the tracer distribution at t = 500 highlights how the material transport resulting from the vortices changes with grass buoyancy for β = 0.06 and β = 0.14 in Fig. 7(a, b).The cores of the vortices lie predominantly above the canopy and the majority of material entrained from the seagrass bed by the vortices appears to come from the upper region of the grass bed (Fig. 7(b)) despite the vortex velocity signature extending well into the grass.Iso-vorticity contours highlight the alignment of the material vortices (Figs 7(a, b)) and the cores of highest vorticity.
Less buoyant blades, β = 0.06, in Fig. 7(a) allow for a larger mean grass deformation and smaller vortices.As the vortices propagate down the channel in Fig. 7(a), they grow in size and in the amount of material entrained.At the same time, material from the overflow is entrained into the grass at a greater rate as well.Vortices are shed from a more or less fixed location x ≈ 15 and grow primarily in vertical extent, as the wavelengths appear to be approximately constant.For β = 0.14 in Fig. 7(b), with more buoyant blades, the vortex size is effectively constant throughout the same domain, meaning that the development region from no vortex to fully developed vortex is shorter compared to the previous case.
To quantify the vertical flux of tracer φ(x, z, t) and the tracer exchange induced by these instabilities, we define the vertical tracer flux as the product of the local vertical velocity and the tracer perturbation with respect to C 0 , expressed as φ = w C , where C = C(x, z, t) − C 0 (z).The instantaneous vertical flux at t = 500 in Figs 7(c, d) for β = 0.06 and β = 0.14, corresponding to the tracer fields in Figs 7(a, b), is positive or negative depending on the co-variance between w and C ; positive flux results from the upward (and downward) movement of anomalously high (and low) tracer anomaly.The flux both out of, and into, the grass grows as the vortex propagates down the channel for smaller β (Fig. 7(c)) .While an individual vortex is experiencing progressively more exchange as it propagates down the channel, the domain as a whole has reached steady-state with regards to the amount of exchange taking place.For β = 0.14 (Fig. 7(d)) the vortices are fully developed and have a constant size and the exchange into and out of the grass bed is balanced.The overturning that occurs inside the vortex cores is reflected by the red-green lobe patterns where the values of φ are dominated by the vertical velocities.
To quantify the relative amount of tracer exchange occurring between the seagrass domain and the overflow, we define the tracer exchange Φ at z = h g , for x ∈ [x a , x b ], as (2) Both φ and Φ should be viewed as relative, as their value depends on the initial tracer distribution.The time evolution of Φ(t) is plotted in Fig. 8(a) for β ∈ [0.02, 0.20], with x a = 20 and x b = 35.The small oscillations observed in each of these curves relates to the vortex turnover time (≈ 10) and to new vortices entering and leaving the domain of integration (see Supplementary information for a video showing how the flux field φ synchronizes with the time-evolution of Φ for β = 0.06 and β = 0.14).Focusing of the long term variations, we observe that the exchange Φ grows once the instability starts, and plateaus in all cases for t > 300.We therefore define the long-term average exchange Φ * as the time-average of Φ(t) for t ∈ [300, 500].
Tracer exchange, measured by this metric, is higher for less deformable blades (larger β), when all other parameters kept constant.Φ * varies with the buoyancy parameter β and with Re (Figs 8(b, c)).The increase in exchange with β and Re eventually saturates for β = 0.10 and Re ≥ 1300.For larger Re, we observed a slight decrease in Φ * alongside an increase in the uncertainty associated with the long-term rms value.
The parametric study previously presented is used to investigate the dependence of the long-term average exchange Φ * on (Re, β).The exchange in Fig. 8(d) shows strong correlation with the w rms trends previously observed in Fig. 6(a).This result is in agreement with the comment in Nepf and Ghisalberti [2008] that the exchange of a scalar would follow the same trend of the exchange of momentum, and decrease as canopy deformability and motion increases.For growing Re or β, however, larger uncertainty is observed.This is also seen in Figs.8(b, c) for large β and Re, respectively, and is potentially related to vortex merger events that become more common with increasing Re and induce temporal fluctuations on the exchange (Supplementary information).

Discussion
In our model, where buoyancy and fluid drag determine the grass blade shape, we find that the Kelvin-Helmholtz-like flow instability weakens with more deflection of the blades.We hypothesize that this result will hold regardless of whether the grass deflection is restored by bending stiffness or buoyancy.As our blades stand up by buoyancy and have no flexural rigidity, the buoyancy parameter β controls the degree of deformation of the grass blades and yields analogous results to the Cauchy number [Luhar and Nepf, 2016;Wong et al., 2020] for flexible beam models (Supplemental information).This hypothesis is further supported by similarities between experiments featuring flexible blades with known elasticity [Ghisalberti and Nepf, 2006] and our model.Changing the deflection mechanism, therefore, should not impact the qualitative behavior.
Schematics from previous studies indicated that clockwise-vortices forward-deflect the grass blades immediately below [Ghisalberti and Nepf, 2002;Nepf and Ghisalberti, 2008;Nepf, 2012;Okamoto et al., 2016;Wong et al., 2020].We find to the contrary that the vortices straighten the grass directly below their core because the horizontal velocity perturbation induced by the vortices acts in the opposite direction to the mean flow.This results in a lower drag force on the grass below the vortex core relative to steady-state, which makes the grass blades more erect.Our schematic of how the vortices induce seagrass motion (Fig. 5(e)) provides a correction to previous schematics in the literature and is consistent with the experimental measurements in Ghisalberti and Nepf [2006] (figures 7 and 8), which report that the smallest horizontal velocity aligns with where the grass blades are most erect and the highest velocities where the grass is most deflected.Our result that more deformable seagrass blades inhibit tracer exchange is consistent with the observations.Nepf and Ghisalberti [2008] find that exchange of momentum is most efficient for rigid canopies and exchange efficiency decreases as the deflection of the canopy increases.They conjecture that the exchange of a scalar should follow a similar trend and decrease as canopy deformability and deflection increase.
Our model of the instability and seagrass also highlights some phenomena within the canopy that warrant further study.For cases of large λ, where the flow speed within the canopy is much smaller than the overflow, the velocity perturbations induced by the vortices cause flow reversal near the seabed, which results in flow separation points that propagate with the vortices and could export sediment from the seabed.Sediment resuspension related to the presence of seagrass canopies has been observed in the field [Adams et al., 2016] and quantified in laboratory experiments, where canopies were found to increase seabed sedimentation compared to bare substrates, and the more blades per unit area, the greater the amount of sediment deposited on the seabed [Barcelona et al., 2021].Another feature that is considered in our model is the free surface and its variation relative to the position to the vortex.While Mandel et al. [2019] have measured experimentally the surface signature of the shear-instability that develops as flow moves through a canopy of rigid rods, their study does not consider the relative phase of the surface signature and the induced vortices.However, they provide a schematic indicating the wave crests immediately above the vortices, which is 180°out of phase with our model.Further experimental investigation of the relative phase as well as an analysis of the free-surface signature of monami with a moving grass bed could provide insight into the potential of remotely observing monami.
Despite qualitative agreement with experiments that feature larger scale oscillations of the grass, our model produces small amplitude blade oscillations.Higher order effects that have been neglected in our model, such as grass inertia, added mass, and virtual buoyancy [Luhar and Nepf, 2016;Wong et al., 2020], can be incorporated to more accurately model the grass meadow for oscillations of higher amplitude.Inertia may introduce another characteristic frequency to the oscillatory motion, and tracking how the blade moves in time and using the instant relative velocity between fluid and blade would more accurately represent the drag for faster, higher amplitude grass motion.Further refinement of the method used to distribute the blade forces onto the computational grid to account for the position of the blade and the center of neighboring cells in both x and z direction is also desirable, as it allows for a more accurate description for higher amplitude oscillations.
While we focused exclusively on the effects of Re and β in this study, variations of the dimensionless parameters r, λ, and F r also impact the results and can be addressed using the current version of the model.Additionally, the model can handle variations in the canopy to free-surface height ratio r and spatially uneven λ(x, z).It can be used to study variable blade number per unit area N , spatially variable grass parameters b, d, and drag coefficient c D .In the current study we explored the two-dimensional vortex regime, but performing three-dimensional simulations with the model would be especially beneficial to study instabilities along the y-direction and vortex interactions.Our model could also be used to compare with field studies of seagrass meadows that quantify fluid exchange above and within the canopy [Hansen and Reidenbach, 2017].Additionally, a study of tracer and sediment transport could be undertaken from a Lagrangian perspective, by applying coherent structure detection methods.This could be used to more conclusively explain how the grass motion impacts material exchange.Finally, applications are not constrained to aquatic vegetation, and our model can simulate atmospheric flows through forests by adjusting the dimensionless parameters to produce the typical canopy deformations and velocity magnitudes.

Conclusions
Our two-phase model of buoyant, deformable, non-shear resistant seagrass blades captures the interaction of flow and submerged canopies, yields shear-instabilities that evolve into vortices and induce an oscillatory motion of the grass blades.While previous schematics of the vortex-grass interaction feature the greatest deflection of the grass immediately below the vortices, our model demonstrates that the velocity perturbation induced by these clockwise vortices acts to make the grass immediately below the vortex more erect than the surrounding canopy, forming the maxima in canopy height.Perturbations induced by the vortices to the background flow increase deflection ahead and behind the vortex.
A stability study of the system as a function of (Re, β) demonstrates the onset of instability.The vertical velocity induced by the instability increases with both Re and β.As Re increases, the shear layer strength increases resulting in stronger vortices.Increased β, corresponding to more buoyant grass that is less deformable, also produces stronger vortices, indicating that the deformability of the grass reduces the vortex strength and delays the instability onset.A scalar field advected with the flow is used to quantify material exchange between seagrass and overflow.Tracer exchange is a function of Re and β, with less deformable blades or larger Re leading to an increased shear above the canopy that results in stronger vortices inducing more exchange.Grass deformation, therefore, inhibits fluid exchange by decreasing the shear magnitude at the canopy top, and therefore the resulting vortex sizes and induced vertical velocities.

Non-hydrostatic model formulation
Numerical simulations are run using a version of the Process Study Ocean Model (PSOM) [Mahadevan et al., 1996a,b], a finite-volume, non-hydrostatic model that we modify to account for the seagrass drag, recompute the blade shapes at each time step as a function of the velocity field, and solve the two-way coupled system of equations in dimensionless form.The governing equations for the fluid are the incompressible Navier-Stokes equations with an added body force term that models the seagrass drag on the fluid and is exclusively applied within the seagrass phase.Following the process of homogenization in Wong et al. [2020], this term accounts for the effects on the flow from all resulting forces applied by multiple seagrass blades on the fluid.The free-surface h satisfies the integral form of the kinematic condition, and a scalar tracer field of concentration C is advected with the flow.
We consider a dimensionless system of equations similar to the one used by Wong et al. [2020], with a critical difference being the inclusion of the buoyancy parameter β to restore the grass deflection instead of the Cauchy number for bending stiffness.The five dimensionless parameters Re, β, r, λ, and F r defined in (1) uniquely determine the flow characteristics and can be varied independently by tuning the dimensional parameters ν * , d, , and N , respectively.The resulting dimensionless governing equations, obtained using the characteristic scales For all the simulations presented, the steady-state velocity profile u(z) is imposed at both the inlet and outlet, and a restoring term is applied to a buffer region at the outflow boundary (omitted in Fig. 1) for x ∈ [L, L r ].The buffer is used to suppress the vortices so the flow matches the outflow boundary condition, thereby minimizing reflections from the boundary.The bottom boundary has a no-slip condition for the velocity (u = w = 0).Pressure is constant at the free-surface.In this study, we explore two-dimensional solutions by allowing no variations in the y-direction within the model.

Buoyant blade equations
The seagrass bed is modeled using a single seagrass blade representative per cell center used in the numerical simulations (with the i-th representative rooted at (x = x i , z = 0), Fig. 9(a)).There are N grass blades per unit area, and each representative models the local averaged blade shape and contribution to the flow.Neglecting flow in the y direction, the motion of each blade representative is confined to the the xz plane, and blade representatives are uniformly distributed in the x-direction.
All blade representatives are inextensible and have constant length .We solve for their shape as a function of the grass buoyancy and fluid load due to drag.The shape of a blade is described by the coordinates x g (s) = (x g (s), z g (s)), which are measured with respect to its base at x = x i and parameterized by the distance s along the blade (s = 0 at the base and s = at the tip, Fig. 9 At every instant t, we assume the grass blades are in equilibrium with the flow, which corresponds to neglecting the blade inertia (note that typically d b and the blade acceleration is small, which makes the inertial term negligible compared to the drag and tension contributions).The blades are buoyant (ρ g < ρ) and do not resist shear (their flexural rigidity EI is negligible) Under these assumptions, the only three forces acting on the blade are: tension, drag, and buoyancy (Fig. 9(b)).Other force terms affecting the blade motion [Luhar and Nepf, 2016], such as virtual buoyancy and added mass, are neglected under our assumptions.The tension T = (T x , T z ) is oriented along the blade.We assume that drag acting tangential to the blade is negligible.The drag per unit length f D is normal to the blade and obeys a quadratic drag law, and the buoyancy per unit length f B points upward: where n = (− cos θ, sin θ) is the upstream normal vector to the blade and z = (0, 1) is the unit vector pointing upward.
After non-dimensionalizing the tension, drag, and buoyancy, the force balance for the blade element becomes The boundary condition T(s = ) = 0 at the grass tip allows us to solve for T(s) by integrating (5) from tip to base.Because T • n = 0, the local blade angle is under the assumption that the blade does not overturn (T z > 0 and |θ| < π/2 along the blade).Finally, integrating dx g ds = sin θ, dz g ds = cos θ (7) from the root up (for s = 0 to ) uniquely determines the instantaneous blade shape if the velocity field is known.

Fluid-blade coupling
The two-way coupling in our model accounts for the impact of the fluid velocity on the shape of the grass, as well as that of the grass shape on the fluid, via the drag force.The canopy height h g (x, t) separating the two phases -seagrass and overflow -depends on the instantaneous positions of all blades.Once the coordinates for all blade tips at the instant t have been determined, h g (x, t) is obtained by fitting a spline through the blade tips.
The relationship between the drag force per unit length f , that acts on a blade representative and is used to solve for the blade shape, and the drag force per unit volume F, that acts on any point in the fluid and appears in the momentum equations, is obtained through the process of homogenization presented in Wong et al. [2020].The force F on the fluid is proportional to λ and to the secant of the local blade representative angle θ, which physically accounts for an increase in the effective number of blades per unit area when neighboring plants tilt.Everywhere within the grass bed, where z ≤ h g (x, t), F = −λ sec θ f , and at the overflow, where z > h g (x, t), we set F = 0 .
Note that while F is defined at any location (x, z), f and θ are evaluated along the blade representatives only and are a function of (i, s).In order to distribute the drag from each blade element to its neighboring cells, we use a Gaussian kernel in the x-direction, with standard deviation 0.3∆x.The relationship between f and F couples the equations for fluid and grass blades, and an iterative under-relaxation method is used at each time step for each grass representative to attain convergence to the equilibrated grass shape.

Numerical grid and conformal map
A conformal map for the vertical grid coordinate accounts for variations in time of the seagrass height h g (x, t) and free-surface height h(x, t), allowing us to reproduce the monami dynamics in an open channel while maintaining a uniform grid in the transformed space.This requires transforming the equations and boundary conditions to solve them in the computational domain, as time and the physical domain evolve [Mahadevan et al., 1996a,b].
We discretize the physical domain with a smooth boundary fitted curvilinear grid, and map this domain onto a computational grid that is rectangular, uniform, and has N i by N k grid intervals in the x and z directions, respectively.The free-surface (z = h(x, t)) is mapped onto the top boundary of the rectangle in the computational domain.The seabed (z = 0) is mapped onto the bottom boundary of the rectangle.Finally, the top of the seagrass bed (z = h g (x, t)) is mapped to the top edge of the N kg -th cell row in the computational domain, so that the bottom N kg cell layers correspond to the seagrass phase, where drag is applied, and the top (N k − N kg ) cells correspond to the overflow phase, where there is no drag.Cells corresponding to each of the two phases are illustrated in green and blue in Fig. 1.
The simulations presented here use N i = 216, N k = 48, and N kg = 24.The grid spacing in the physical domain is ∆x = 0.2 horizontally, so that L r = 43.2(the physical channel is 43.2 by 1), and ∆z ≈ 0.02 vertically (note that ∆z is non-uniform and varies depending on h g and h), and a uniform time step ∆t = 0.1 is used.The horizontal length of the domain before the buffer region where the velocity profile is restored to u(z) is L = 36, and the dimensionless domain considered for the analysis is

A Comparison between buoyancy and rigidity models
Prior to the onset of instability, or closer to the inflow section of the channel, the flow and grass are both steady.The steady state solution is a function of z alone, and can be calculated with a simplified onedimensional coupled model that eliminates dependence in x and t (Methods Section).
The steady state velocity profiles u(z) and corresponding grass positions (x g , z g ) calculated for Re = 10 3 , r = 0.5, λ = 1, F r 2 0.1, are presented for a range of values of the buoyancy parameter β (Fig. S1).Solutions are computed for two different fluid boundary conditions: du/dz = 0 (free-slip) at the bottom (z = 0) (Fig. S1(b)) and u = 0 (no-slip) (Fig. S1(c  Grass blades with bending stiffness (Fig. S1(a)) and buoyancy (Fig. S1(b)) result in qualitatively similar solutions, with the main difference arising from the clamped bottom boundary condition in Wong et al. [2020] that prevents deflection, compared to the hinged bottom boundary condition we used.The hinged condition results in lower drag near the bottom boundary.Another difference is the fluid velocity at the tip of the grass.For the buoyant model in Fig. S1(b), we observe a monotonic growth of the velocity at the tip with decreasing β, while in Fig. S1(a) there is a peak and then a decay for increasing C Y .Additionally, small variations of β, as β → 0, drastically change the steady-state solutions (note the differences for β = 0.01 and 0.015 in Fig. S1(b), for example).
The no-slip bottom boundary condition (Fig. S1(c)) reduces blade deflection at the root, resulting in a smaller range of grass deflection angles and ultimately less deflection at the tip.The shape of the velocity profile near the canopy top and the velocity shear magnitude are not sensitive to the bottom boundary condition choice.The no-slip condition at z = 0 is used for the time-dependent simulations as it is physically more accurate.
B Vortex merger events and tracer exchange At combinations of high Re and high β, larger-scale time fluctuations in Φ(t) become more apparent (Fig. 8).
To understand what contributes to these additional time scales to the flow, we study the vortex structure and tracer concentration fields for Re = 1500.The simulations presented in Fig. S2 only vary the buoyancy parameter, with the left and right panels corresponding to β = 0.06 and β = 0.20, respectively.Figs S2(a, b) present the vorticity ζ, and Figs S2(c, d) the tracer concentration C, both at time t = 300.These plots highlight how the vorticity field lines up accurately with the tracer field at a given time, even though the tracer concentration has been evolving for t ∈ [0, 300].
While for β = 0.06 (Figs S2(a, c)) all vortices look similar and periodic, from both vorticity and tracer perspectives, the β = 0.20 results (Figs S2(b,d)) show signals of vortex interaction, with an imminent vortex merger [Ho and Huerre, 1984] that has started at x ≈ 26 (see Supplementary information for a video of the time evolution of these merger events).
Figs S2(e, f ) present Hovmöller diagrams of w(x, z = h g , t) for both cases, similar to the ones in Figs 4(a, b).A periodic signal with vortices propagating at constant speed and amplitude is observed for the more deformable blade case (Fig. S2(e)), while more nonlinear interactions and vortex merger events are observed for the less deformable case (Fig. S2(f )).Both cases have a similar dominant vortex speed that is again close to the 0.6 obtained for β = 0.10 in Fig. 4(b).However, the w amplitudes for β = 0.20 are stronger, more variable, and the periodicity of the signal is less apparent.This variability may be due to bigger vortices for β = 0.20 growing enough to interact with neighboring vortices, causing merger events and other nonlinear phenomena that reduce periodicity in the flow.No preferred frequency for vortex merger events was observed.

C Video description
The four supplementary videos, for which captions are available below, correspond to the time evolution of the results presented in Figs 3, 7, 8, and S2.
Video 1 Instability onset.(top) Vorticity field ζ for the full domain.(bottom) From left to right, steady-state horizontal velocity u, horizontal velocity perturbation u , vertical velocity perturbation w, and grass blade positions (x g , z g ), for the designated region of the domain (dashed rectangle).Black lines represent the instantaneous seagrass height h g (x, t).Video corresponds to Fig. 3(a, b) in the manuscript.

Figure 4 :
Figure 4: Hovmöler diagrams presenting the space and time evolution of the vertical velocity at the top of the grass bed w(x, z = h g , t) for (a) the full domain and (b) a focused region indicated by the black box in (a).The yellow line in (b) indicates a constant speed 0.6 of propagation of the perturbations.(c) Vertical velocity rms along x at t = 0, 50, and 500.(d) Steady-state background horizontal velocity.(e) Time evolution of the spatial rms of w (blue) and (h g − h g ) (green), identifying the instability onset and long-term behavior from both fluid and grass perspectives.

Figure 5 :
Figure 5: Instantaneous plots at t = 500, for the (a) surface height h (solid blue line) and grass height h g (solid black line, with dashed line representing the reference steady-state height h g ), (b) contours of vorticity ζ and arrows representing the velocity field u, with the grass height perturbation (black line) amplified by a factor of 4. (c) Perturbation velocity field u = u−u and grass blade representatives, with colors representing the angle of deflection with respect to the steady-state angle, δθ = θ − θ, in radians.(d) Histogram of the grass height distribution within the entire domain, with inset representing the blade shapes corresponding to minimum (red), maximum (blue), and steady-state (black) height (bin width = 0.001).(e) Schematic of how vortices (iso-vorticity contours) induce grass blade deflection, with arrows representing the velocity perturbation induced, and colors indicating how they deflect the grass blade locally.The blue line on top representing the free-surface.

Figure 6 :
Figure 6: Parametric sweep of (Re, β) combinations and the resulting long-term rms of the vertical velocity and grass height perturbation, for r = 0.5, λ = 1, and F r 2 = 0.1.Contour plots for (a) w rms tracking instability from the flow perspective, and (b) δh rms g from the seagrass bed perspective.Lines represent contour levels.

Figure 8 :
Figure 8: (a) Tracer exchange Φ as a function of time for β ∈ {0.02, 0.04, . . ., 0.20}.Long-term exchange Φ * (b) as a function of β, for Re = 1000, (c) as a function of Re, for β = 0.10, and (d) as a function of (Re, β).Error bars in (b, c) correspond to the rms deviation from the mean, and colors match curves in (a).In (d), the colors represent the mean value Φ * for the given parameter combination, and the marker size represents the rms deviation from the mean.
(a)).The blade coordinates in the xz plane are uniquely determined by x i and the clockwise blade angle θ(s) with the vertical.The blade has width b (along y) and thickness d (along x).

Figure 9 :
Figure 9: (a) Schematic of a grass blade representative, rooted at x i , with grass coordinates in green.(b) Forces acting on an infinitesimal, inextensible blade element: tension T, drag per unit length f D , and buoyancy per unit length f B .
)).At the surface (z = 1), du/dz = 0.As a comparison, Fig.S1(a) presents the solution ofWong et al. [2020] for neutrally-buoyant blades with flexural rigidity EI and a free-slip bottom boundary condition, where the parameter controlling the blade deformability is the Cauchy number C Y = ρbC D H 3 U 2 /EI.In our model, EI → 0, C Y → ∞, buoyancy is the dominant agent that resists the fluid drag on the blade, and therefore β quantifies to what extent the blade can deform.For these steady-state solutions, β is varied while keeping Re, λ and r constant.Note that β and Re can be adjusted without modifying any of the other dimensionless groups by tuning the blade dimension d and ν * , respectively.

Figure S2 :
Figure S2: Vortex interaction for variable β and fixed Re = 1500.Instantaneous (a, b) vorticity ζ(x, z, t = 300) and (c, d) scalar field C(x, z, t = 300).(e, f ) Hovmöller diagrams presenting the space and time evolution of w(x, z = h g , t) for β = 0.06 ((a, c, e), more deformable) and β = 0.20 ((b, d, f ), less deformable), highlighting vortex interactions and merger events for β = 0.20.The yellow slope indicates the reference dx/dt = 0.6, and the dashed black line marks the time t = 300 when the fields are plotted in (a − d).