Wind exposure and sediment type determine the resilience and response of seagrass meadows to climate change

Seagrasses and bare sediment represent alternative stable states, with sediment resuspension being a key driver of system stability via the Seagrass–Sediment–Light (SSL) feedback. We explore the SSL feedback by quantifying the sediment stabilization by seagrass, and using these measurements to calculate under which conditions seagrass ends up in a turbid environment. We quantified in‐situ sediment resuspension velocity thresholds (ucr) for Zostera marina growing in medium to fine sand, using a field flume inducing near‐bed wave motion. ucr was determined for full length shoots, shoots clipped to 0.08 m, and removed shoots. We found that rhizomes did not influence ucr of the top sediment layer. Overall, ucr was linearly related to blade area, which became independent for sediment type when normalizing ucr for the resuspension threshold after shoot removal. Comparing measured ucr against natural wave conditions showed that the seagrass meadow at the study site is currently stable. Exploring the effects of changing hydrodynamic conditions revealed that effects of increasing storminess has limited influence on sediment resuspension and thus the SSL‐feedback. Increasing mean wind velocity had a stronger influence on SSL‐feedback dynamics by causing more frequent exceedance of ucr. The response of seagrasses to increasing wind pressure depends on bay topography. A fully exposed Z. marina meadow under low initial turbidity pressure trended toward bistability, as turbidity pressure increased mainly on bare sediments. The study site and a fully exposed Z. marina meadow under high initial turbidity pressure saw an increase in turbidity across all blade areas.

Seagrass meadows are among the most productive and valuable marine ecosystems (Costanza et al. 1997;Duarte and Chiscano 1999). The structural complexity of seagrass beds makes them suitable as a nursery (Heck et al. 2003, and references therein), and they are important marine carbon sinks (Duarte and Cebri an 1996). Seagrass meadows also benefit coastal ecosystem stability by their ability to attenuate hydrodynamic energy from both currents and waves (Gambi et al. 1990;Fonseca and Cahalan 1992;Infantes et al. 2012). Locally this reduction in hydrodynamic energy allows seagrasses to trap sediment and prevent erosion during stormy conditions as rhizomes fixate the sediment trapped by the blades, reducing sediment erodibility (Christianen et al. 2013). This makes seagrass meadows beneficial for maintaining shorelines by reducing beach erosion (James et al. 2019). Because of these local and remote sediment stabilizing properties seagrasses are classified as ecosystem engineers, that can be used for nature based erosion and flood protection (Ondiviela et al. 2014).
Seagrass canopies reduce hydrodynamic energy, facilitating sediment trapping by reducing resuspension under stormy conditions (Gacia and Duarte 2001). This can decrease the water turbidity locally, which results in a higher number of plants as higher light availability favors growth. These can in turn retain more sediment. Overall, this results in a positive feedback loop. From now on, we refer to this as the seagrasssediment stabilization-light (SSL) feedback.
Positive feedback loops are also known to induce bistability dynamics (Scheffer et al. 2001;van der Heide et al. 2007). In *Correspondence: jaco.de.smit@nioz.nl This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Additional Supporting Information may be found in the online version of this article.
Special Issue: Nonlinear dynamics, resilience, and regime shifts in aquatic communities and ecosystems seagrass meadows, this may be manifested by a period of increased pressure causing biomass loss, e.g., due to hydrodynamic or anthropogenic disturbances, which may lead to increasing turbidity and the subsequent collapse of seagrass meadows and thereby the SSL-feedback (Yaakub et al. 2014;Maxwell et al. 2017;Adams et al. 2020). To recover, the hydrodynamic energy needs to reduce to a lower level than when seagrass was present, in order to have sufficiently low turbidity levels that enable seagrass establishment in the absence of the SSL-feedback (van der Heide et al. 2007). Bistability and threshold behavior leads to low success rates of seagrass restoration efforts, as successful restoration often requires transplantation of large numbers of shoots or bio-mimicry to overcome establishment thresholds (van Katwijk et al. 2016;Temmink et al. 2020). The nonlinear responses that are typical for alternative-stable state dynamics make it difficult to predict at which point a seagrass meadow will collapse during increasing hydrodynamic disturbance.
Managing seagrass beds for nature values, erosion protection, or as part of a flood defense scheme thus requires quantitative understanding of both sediment erodibility in seagrass meadows, and seagrass persistence under threshold hydrodynamic conditions. Global change may affect the prevalence of these conditions by changing sea level and wind conditions, which are projected to cause wave climate change in approximately 50% of the global coastline (Morim et al. 2019). Threshold conditions are exactly the conditions where seagrasses are most needed for the valuable services they provide, but also where they may be at the highest risk of collapse. Despite recent studies on the physical mechanisms determining sediment stability in vegetated canopies (e.g., Hansen and Reidenbach 2012;Yang et al. 2016;Tinoco and Coco 2018), there is a significant knowledge gap on (especially waveinduced) resuspension thresholds in natural seagrass meadows, and how these are affected by the physical properties of the seagrass and the grain size distribution of the sediment (De Boer 2007). Maximum unidirectional flow thresholds for seagrass presence have been determined for some species, and vary between 7 and 150 cm s À1 (Koch 2001, and references therein). There is however no clear distinction in the literature between mechanical thresholds such as uprooting and cliff formation by edge scouring, and sediment resuspension thresholds disturbing the SSL feedback. This lack of quantitative data is due to the fact that field measurements do not allow for control of wave properties, and laboratory flume studies typically lack natural meadows on undisturbed sediment.
In this study, we therefore used a wave-generating field flume (the TiDyWAVE flume; de Smit et al. 2020) to (1) measure in-situ resuspension thresholds inside a Zostera marina meadow and (2) to quantify the dependence of such resuspension thresholds on the combined effects of shoot length-density combinations, presence of rhizomes, and sediment type. Subsequently, we (3) use a simple wave generation model to compare our measurements with natural wave conditions and (4) to assess to what extent altered hydrodynamic pressure due to changing wind conditions affect turbidity pressure.

Study-site characteristics
Resuspension threshold measurements were conducted in Bokevik bay, located in the Gullmars fjord next to the Kristineberg Research Station, Sweden (Fig. 1a,b). The bay is characterized by sheltered conditions during the dominant south-westerly winds. During north-easterly winds, the maximum fetch is 19 km long, leading to wave heights of approximately 0.5 m with peak orbital velocities in the order of 0.55 m s À1 inside the seagrass meadow during storms. Inside the bay, the average water depth is 0.5 m, with a 0.1 m tidal amplitude. Z. marina occurs in a narrow band at the seaward edge of the shallow bay, between 1 and 4 m depth (Fig. 1c).
Two locations at the shallow edge of the seagrass meadow were selected for the resuspension threshold measurements based on their contrasting seagrass morphology and sediment type ( Fig. 1c and Table 1). Site LoFi (i.e., Long shoots and Fine sediment) is characterized by longer Z. marina shoots (0.24 AE 0.04 m canopy height), which was measured as the distance from the bed to the top of the canopy. The sediment, which was analyzed with a Malvern laser particle sizer, had a mean grain size of 135.51 AE 4.21 μm and a 5.26 AE 2.09% mud content. The shoot density was 262 AE 38 shoots m À2 , which was measured by counting the number of shoots inside the TiDyWAVE flume (i.e., over a 0.6 Â 0.4 m area). Site ShoCo (i.e., Short shoots and Coarse sediment) is characterized by shorter shoots (0.17 AE 0.02 m canopy height). The sediment was coarser than LoFi with a mean grain size of 244.24 AE 6.90 μm, and 0% mud content. The shoot density was however similar with 295 AE 48 shoots m À2 .
The blade area per square meter of seabed area, which causes the drag on water flow and thus the protection of the sediment underneath, could not directly be determined within the flume experiment plots because blades were clipped and removed during the field flume experiments. Therefore, four biomass cores (0.25 m diameter) were taken at each site. The blade lengths, blade width, number of shoots, and canopy height were measured of the seagrass from the biomass cores, yielding individual blade areacanopy height Â shoot density relations for site LoFi and ShoCo (Fig. S1). These were used to calculate the total blade area in the flume experiment plots using the measured shoot density and canopy height.  Table 1). near-bed flow of natural waves in terms of velocity and turbulence characteristics (de Smit et al. 2020). The wave period was set at T = 4 s, which corresponds to typical wave periods during stormy conditions in fetch-limited bays such as Bokevik bay. The TiDyWAVE flume was placed carefully over the seagrass bed, and oscillatory flow was gradually increased until incipient motion was observed visually by trained divers. Incipient motion was defined as the frequent movement of sediment particles over the bed across the entire measurement section. The incipient motion threshold was interpreted as the resuspension threshold, as when incipient motion and ripple formation occurs the finer sediment particles and lower-density organic matter present in the sediment will be suspended in the water column. This will consequently lead to light attenuation, which will affect the SSL-feedback. Effects of, e.g., bioturbation or sediment binding by benthic microalgae are thus not accounted for explicitly, but are included implicitly through conducting these measurements in situ.
To determine the effect of shoot length and the effect of rhizomes only, resuspension thresholds were determined for (1) fulllength shoots, i.e., 0.24 and 0.17 m long shoots for respectively LoFi and ShoCo, (2) shoots clipped to 0.08 m length above the sediment, that is, representing Z. marina seedlings (c.f. Orth and Moore 1983), and also obtaining equal shoot lengths for ShoCo and LoFi, and (3) shoots cut at the sediment surface but leaving the rhizome structure intact. These treatments were made at the exact same spot of the initial resuspension threshold measurement to ensure that sediment properties and shoot density remained constant between the different treatments. Sediment disturbance between the experiments was minimized by careful clipping and stopping the experiment within 3 min after incipient motion was observed. Statistical significance of the differences between the shoot length treatments were assessed with paired T-tests, as length treatments were conducted at the exact same spot. Statistical significance of the effect of rhizomes on sediment erodibility was assessed by comparing the resuspension threshold of fully removed shoots and bare sediment outside the seagrass meadow using Welch's T-tests, as these measurements were conducted at different spots.

Hydrodynamic measurements in the field flume experiment
Orbital velocities were measured with an Acoustic Doppler Velocimeter (ADV, Nortek ® Vectrino Profiler) which was inserted into the flume from the side at 0.15 m above the bed. This height was chosen so that the ADV measured the freestream velocity above the seagrass canopy as it bends down during peak velocities. In contrast to near-bed velocity under flexible vegetation canopies, which provides insight in hydrodynamic processes within the canopy, the free stream velocity can be used to calculate the corresponding wave conditions. Hence, the results of hydrodynamic processes within the Z. marina canopy are measured instead of the processes themselves. Given that the canopy of flexible vegetation continuously changes in vertical position, the size of wave boundary layer at the canopy-water interface remains in the order of millimeters and its location varies over time. Hence, the canopy height does not influence the time averaged ADV measurements as long as they are conducted above the deflected canopy height, as there is no discernible boundary layer effect on time averaged measurements (van Veelen et al. 2020). Because peak orbital velocity in TiDyWAVE varied slightly between individual waves, the velocity at which incipient motion was observed, hereafter referred to as u cr , is defined as the mean of all oscillatory flow velocity peaks during the period that incipient motion was observed. To remove disturbances in the ADV data due to interference from, e.g., seagrass blades and algal matter, data were filtered for signal to noise ratios above 30, signal amplitude values above À40 dB, and pulse-reflection correlation values above 90%.
Calculating exceedance probability values of measured u cr under present and potential future wave conditions

Calculating near-bed orbital velocity
Local exceedance probabilities of u cr were calculated using the Jonswap method. While this method does not solve important 3D processes like depth dependent light attenuation and presence of suspended sediments from nonlocal sources, it does provide a basic indication of present turbidity pressure inside the Z. marina meadow and allows for conducting a sensitivity analysis to changing wind conditions. The Jonswap method uses wind velocity, duration, and fetch to calculate the corresponding wave heights and periods. The underlying assumption of this method, that wave generation is not limited by depth, is valid as the depth in the center of the Gullmars fjord ranges from 45 to 125 m. The Jonswap method can also be used to calculate the amount of time it takes to generate a fully developed sea where waves are fetchlimited. In the case of Bokevik bay it takes less than 3 h to obtain a fully developed sea for wind velocities above 13.9 m s À1 (i.e., stormy conditions) and a fetch below 19 km. So, for simplicity, waves were assumed to always be fetch-limited. The wave height and wave period can then be calculated as: where H is wave height (m), T is wave period (s), U is wind velocity at 10 m above ground level (m s À1 ), g is gravitational acceleration (m s À2 ) and H* and T* are dimensionless wave height and period, calculated as: de Smit et al. Seagrass resilience to changing winds S124 where F * is the nondimensional fetch, which is calculated as: where F is the fetch (m). From the wave height the corresponding near-bed peak orbital velocities can be calculated using linear wave theory, which is not strictly valid for shallow water, but gives reasonable estimates when applied in the shoaling zone where the maximum error is approximately 20% (Guza and Thornton 1980). Peak orbital velocity is calculated as: where u is the peak near bed velocity (m s À1 ), d is water depth (m), ω ¼ 2π T is the angular velocity of the wave (rad s À1 ), and k ¼ 2π L is the wave number (m À1 ) where L is wavelength (m). The wavelength was iteratively calculated from the dispersion relation: Model calibration The Jonswap model was calibrated using ADV measurements collected from 3 to 17 November 2016, and corresponding wind data downloaded from the nearest wind gauge of the Swedish Meteorological and Hydrological Institute, located 18.5 km from the study site (Fig. 3b). The ADV measurements were taken inside the seagrass meadow at 0.15 m above the seabed, similar to those of the field flume experiments. The wind gauge measured at 15 m above ground level, so wind velocity was corrected to 10 m above ground level assuming a logarithmic velocity profile: where U 10 and U 15 are wind velocity at respectively 10 m and 15 m above ground level. The mean water depth was 2 m. During the measurement period both mild and strong winds occurred from both the southwest and northeast. These are the conditions of most interest as NE winds generate the strongest waves, and SW winds are the dominant winds in this area. The narrow, elongated shape of the Gullmars fjord was observed to cause unrealistically strong predicted changes in wave height with minor changes in wind direction. Therefore the wind direction -fetch relation was smoothed based on back-calculated fetch from the ADV data (Fig. S2, Fig. 3a,b). As the Jonswap model performed poorly under calm conditions (Fig. 3c), i.e., well below the resuspension threshold of the bare sediment, peak orbital velocities below 0.1 m s À1 were not considered in the assessment of the calibration. Because of the simplified assumptions in the model the root mean square error of the calculated orbital velocities was 8.1 cm s À1 , increasing toward lower orbital velocities (Fig. 3d). However, the model is able to predict the general distribution of orbital velocities and is more accurate toward higher orbital velocities, i.e., the hydrodynamic pressure to which the Z. marina meadow is exposed (Fig. 3d).

Appending climate change scenarios
The influence of changing hydrodynamic conditions due to climate change on the stability of Z. marina in Bokevik bay was assessed by changing the wind conditions in the wave model, and calculating the resulting exceedance probability, P (u > u cr ), of the measured u cr . In other words, the probability that water is turbid for these future scenarios over a range of Z. marina blade areas was calculated. Sea level change was not considered, as this generally results in lateral migration of seagrass meadows rather than collapse (Duarte 2002). Wind scenarios were applied to Bokevik bay and to two fully exposed conceptual bays, where fetch is not influenced by wind direction. The conceptual bays consisted of a bay with a 1.9 km constant fetch (C 1.9 ), yielding an initial turbidity pressure inside the full canopy Z. marina meadow similar to that of Bokevik bay, and a bay with a 6 km constant fetch (C 6 ), yielding a high initial turbidity pressure, near the collapse threshold of approximately 50% light limitation (Biber et al. 2009), inside the full canopy Z. marina meadow. This allows for (1) assessing the effect of wind sheltering due to bay topography, and (2) assessing to which extent wind climate change effects differ along a turbidity pressure gradient.
Two wind scenarios were considered: (1) increasing mean wind velocity, and (2) increasing storminess, which is defined here as the percentage of winds exceeding 13.9 m s À1 , that is, 7 bft or higher. The wind scenarios were incorporated as an alteration of the present distribution of wind velocities at Bokevik bay, which was derived from wind data of the same weather station as was used for model calibration. In order to account for interannual variability, 8 yr of wind data, ranging from January 2011 until December 2018, was used. Wind velocities follow a Weibull distribution, of which the probability density function reads: The shape parameter a and the scale parameter b can be derived from the wind velocity timeseries. Based on the distribution of real wind velocity at Bokevik bay, a random wind velocity distribution was generated (Fig. S3a). Similar random wind velocity distributions were generated for the increasing mean wind velocity and increasing storminess scenarios by iteratively changing the a and b parameters of the baseline Weibull distribution to yield either a 5%, 10% and 15% increase in mean wind velocity without de Smit et al. Seagrass resilience to changing winds changing storminess, or a 5%, 10% and 15% increase in storminess without changing the mean wind velocity (Fig. S3b).
For the two conceptual bays, the entire wind velocity timeseries was used to derive the wind scenarios. Hence, the resulting wave heights, near-bed orbital velocities, and exceedance probabilities of u cr can be calculated directly. For Bokevik bay however, given that wind velocity and direction are correlated and fetch is variable (Fig. 3a), wind velocity distributions were calculated for 10 wind direction bins, to which the wind climate scenarios were applied individually (Fig. S3c). Wave height, nearbed orbital velocities, and P(u > u cr ) were calculated for each wind direction bin, and multiplied with the occurrence frequency of that wind direction to yield the total distribution:

Resuspension threshold measurements
The flume measurements clearly demonstrate that the presence of Z. marina increased the resuspension threshold velocity (u cr ) at both locations, and that the resuspension threshold is linearly related to Z. marina blade area per m 2 seabed area (Fig. 4a). On ShoCo, Z. marina increases u cr with 26% from 0.21 AE 0.02 m s À1 (mean AE SD) on bare sediment to 0.26 AE 0.01 m s À1 inside the full length meadow. On LoFi, Z. marina increases u cr with 47% from 0.18 AE 0.01 ms À1 on bare sediment to 0.26 AE 0.01 m s À1 inside the full length meadow.
On ShoCo, u cr does not differ significantly (Welch's T-test p = 0.21) between bare sediment adjacent to the seagrass meadow (u cr = 0.21 AE 0.02 m s À1 ) and seagrass where the aboveground biomass was removed (u cr = 0.20 AE 0.01 m s À1 ). This indicates that belowground biomass does not influence u cr of the top sediment layer at this site. In contrast to ShoCo, at LoFi a larger amount of fine sediment is trapped inside the Z. marina meadow compared to the adjacent bare sediment (Table 1). As a result, at LoFi the value of u cr is lower (Welch's T-test p = 0.01) when seagrass was removed (u cr = 0.14 AE 0.01 m s À1 ) compared to the control measurements on bare sediment (u cr = 0.18 AE 0.01 m s À1 ), which did not contain this finer material.

Predicting reduction of u cr and bed shear stress attenuation based on seagrass morphology
Interestingly, even though sediment type and shoot morphology are distinctly different between ShoCo and LoFi, u cr did not differ significantly between the sites with full length seagrass (Fig. 4 a, Welch's T-test p = 0.67). However, the increase in u cr for a given blade area is smaller, indicating that Z. marina is less efficient in retaining sediments toward finer de Smit et al. Seagrass resilience to changing winds grain sizes. When u cr of the full-length and clipped shoot treatments is standardized with u cr when shoots are removed, both curves in Fig. 4a collapse on the same line (Fig. 4b). This indicates that the difference in u cr between ShoCo and LoFi can be explained by the difference in sediment stability, and that the ability of seagrasses to increase sediment resuspension thresholds is a function of blade area, which is independent of sediment type. Given that the critical bed shear stress of the sediment is constant, the reduction in bed shear stress by the Z. marina canopy at the sediment resuspension threshold can be calculated as the square of the standardized reduction in u cr . This yields an asymptotic relation between Z. marina blade area and bed shear stress attenuation by the canopy at the sediment resuspension threshold (Fig. 4c). Low-density seagrass beds are able to significantly reduce bed shear stress, while additional reduction in bed shear stress toward higher blade densities diminishes.

Effects of changing wind conditions due to climate change
In both fully exposed conceptual bays and in Bokevik bay, the effect of increasing storminess on P(u > u cr ) is minimal ( Fig. 5a-c), indicating that SSL feedback dynamics in seagrass meadows are not sensitive to changing storminess. Only the fully exposed bay with equal turbidity pressure inside the Z. marina meadow as Bokevik bay (C 1.9 ) showed a minor increase in P(u > u cr ), approximately 1% depending on sediment type and Z. marina blade area, toward high blade areas (Fig. 5b). This was similar for both ShoCo and LoFi. The limited effect of increasing storminess on P(u > u cr ) can be explained by the fact that winds below the storm threshold are able to generate sufficiently large waves to cause sediment resuspension (Fig. 6). This is especially the case in sheltered bays such as Bokevik bay, where sediment resuspension events are dominated by favorable wind direction rather than high wind velocity (Figs. 6 and S4). Only in high density, fully exposed Z. marina meadows under a low turbidity pressure, sediment resuspension events are storm dominated. However, given that initial P(u > u cr ) under these conditions is low, the effect of increasing storminess remains small.
Increasing the mean wind velocity had a much larger effect on P(u > u cr ) than increasing storminess, and it's effect is contrasting between the fully exposed conceptual bays and Bokevik bay, and contrasting between ShoCo and LoFi de Smit et al. Seagrass resilience to changing winds ( Fig. 5a-c). In Bokevik bay and the fully exposed conceptual bay with high initial turbidity pressure inside the Z. marina meadow (C 6 ), an increase in mean wind velocity leads to a significant increase in P(u > u cr ) on both bare sediments and across the assessed blade area gradient (Fig. 5a,c). Hence, an increase in mean wind velocity causes an increase in turbidity pressure on both bare sediment and inside seagrass meadows in these systems. Similar to the small effect of storminess on P (u > u cr ) in these cases, the larger effect of increasing mean wind velocity is caused by the majority of winds leading to sediment resuspension being below the storm threshold (Fig. 6). In C 1.9 , an increase in mean wind velocity causes a significant increase in P(u > u cr ) on bare sediment and low blade areas, but no increase toward high blade areas (Fig. 5b), as sediment resuspension events at these blade areas are storm dominated (Fig. 6). Hence, the difference in turbidity pressure between bare sediments and seagrass meadows increases, indicating a tendency toward bistability. While increased storminess yielded similar changes in P (u > u cr ) for ShoCo and LoFi, increasing mean wind velocity showed considerable differences. In Bokevik bay and C 6 , the increase in P(u > u cr ) for LoFi was relatively consistent, while for ShoCo it reduced toward higher blade areas (Fig. 5a,c). This indicates that dense Z. marina meadows on coarse sediment are more resilient to increasing mean wind velocity than Z. marina meadows on finer sediment. Similarly, the increase in P(u > u cr ) on bare sediment and low density Z. marina meadows with increasing mean wind velocity in C 1.9 was lower for ShoCo. This indicates that seagrass meadows on coarse sediments are less susceptible to bistability dynamics under increasing wave pressure.

Discussion
Field flume measurements on resuspension thresholds in a Z. marina meadow revealed a linear relation between blade area per m 2 seabed area and the near-bed peak orbital velocity threshold for sediment resuspension (u cr ). When normalizing u cr with u cr of removed aboveground biomass, the blade areau cr curves for both sediment types collapse. This relation is thus independent of sediment grain size and the presence/ absence of rhizomes. The Seagrass-Sediment-Light (SSL) feedback-related seagrass meadow stability was assessed by combining the field flume measurements with a data-driven wave generation model. As the Z. marina meadow in Bokevik bay grows at a very sheltered location, the meadow is not  Fig. 5. Effects of increasing mean wind and increasing storminess (color gradient indicates 5, 10, 15% increase) on P(u > u cr ) along a range of Z. marina blade areas from bare sediment to full canopy for (a) Bokevik bay, (b) a fully exposed bay with similar turbidity pressure in the Z. marina meadow as Bokevik bay (1.9 km constant fetch), and (c) a fully exposed bay with a Z. marina meadow under high initial turbidity pressure (6 km constant fetch). de Smit et al.
Seagrass resilience to changing winds prone to light limitation due to sediment resuspension under both present and potential future hydrodynamic conditions. Only increasing mean wind velocity, and not increasing storminess, had a significant effect on u cr exceedance probability (P(u > u cr )). Wind direction can play an important role in fjords, as their narrow and elongated topography allows only a narrow range of wind directions where fetch is sufficient to generate strong waves. Therefore, calmer winds with maximum fetch more often caused exceedance of the resuspension threshold rather than storms.

Sediment resuspension in seagrass meadows
The field flume measurements showed that the effect of blade area on u cr is weaker for LoFi. More blade area is required to obtain the same increase in u cr compared to ShoCo. The blade area range where Z. marina is potentially unstable thus increases toward finer sediments. This observation is in line with previous modeling results (Carr et al. 2010), which showed that the range of water depths where seagrasses are bistable is larger with finer sediments as the lower erodibility of these sediments leads to increased sediment resuspension in shallow water. The similarity between u cr when shoots are removed and u cr outside the seagrass meadow shows that u cr for sediment resuspension is determined by aboveground biomass rather than rhizomes. Our sediment resuspension measurements are however limited to the top sediment layer only, as the short duration of the experiments did not allow for substantial erosion. Rhizomes still prevent further erosion by stabilizing sediment directly below the top layer (Marin-Diaz et al. 2020), and thus still provide important coastal protection services even with low aboveground biomass (cf. Christianen et al. 2013). However, as initial resuspension of the top sediment layer already leads to increased turbidity, the influence of rhizome biomass on SSL-feedbacks inducing bistability is likely limited.

Effects of shoot density and blade area on near-bed turbulence and u cr
There is a long line of research aiming to obtain accurate measurements and predictions of u cr within vegetation. This is complicated, as under oscillatory flow a shear layer develops at the canopy-water interface, and flow at the bed is reduced (Lowe et al. 2005). However, in case of a short and sparse canopy, turbulent vortices penetrate further through the canopy and reach the bed, increasing near-bed turbulent kinetic energy and reducing near-bed u cr compared to bare sediments (Yang et al. 2016;Tinoco and Coco 2018). Moreover, shorter shoots are less bendable, which further increases turbulence at the bed (Pujol et al. 2013). Therefore, Tinoco and Coco (2018) suggest that near-bed turbulent kinetic energy (TKE) explains the beginning of sediment motion rather than orbital velocity, as the bed shear stress is a direct function of TKE. The results of Tinoco and Coco (2018) confirm that the critical bed shear stress for sediment resuspension is unaffected by the presence of a canopy, as the critical near-bed TKE did not change with shoot density in their experiments. These measurements are valuable from a mechanistical point of view, but they do not provide direct information on critical wave conditions for sediment resuspension. Quantifying critical wave conditions from near-bed hydrodynamic conditions would require additional quantification of the energy attenuation by the seagrass blades. Therefore, we measured the orbital velocity above the seagrass canopy instead. While inherently accounted for, this method does not provide detailed mechanistical information on the hydrodynamic processes within the canopy. It does however provide direct insight in the wave climates under which sediment resuspension occurs, as these can be calculated from, e.g., the Jonswap equations and linear wave theory.
Unlike laboratory studies showing decreasing near-bed u cr , i.e., the near-bed flow velocity needed to generate a certain TKE reduces with increasing shoot density (Yang et al. 2016;Tinoco and Coco 2018), we did not observe a shoot density effect on shear stress attenuation. These studies however used relatively low shoot densities (25-250 shoots m À2 ), and the mimics were rigid. This range of shoot densities is lower than shoot densities generally occurring in seagrass meadows, e.g., $140-2200 shoots m À2 for Z. marina (Olesen and Sandjensen 1994). In natural seagrass meadows, sparse canopies may indeed enhance turbulence, but a reduction is observed toward dense canopies (Hansen and Reidenbach 2017). Moreover, the results of Yang et al. (2016) and Tinoco and Coco (2018) show a decreasing effect of shoot density on u cr reduction toward higher shoot densities. Therefore, any direct effects of shoot density on u cr will be much smaller than the corresponding effect of changing blade area per m 2 in the case of a natural Z. marina meadow.

Assessing seagrass meadow resilience to climate change
The modeling approach used in this study to assess potential changes in future seagrass meadows stability due to changes in hydrodynamic conditions remains minimalistic, as it is a 1-D approach which assumes that water is either fully turbid or fully clear depending on the exceedance of the sediment resuspension threshold. It therefore does not resolve some important 3-D processes such as water residence time (Adams et al. 2018), sediment diffusion and depth dependent light attenuation (Carr et al. 2010), suspended sediments from nonlocal sources (Carr et al. 2018), and bed level change in case of long term predictions. While the current model is able to predict turbidity events reasonably well in a small bay with a local source for suspended sediments, the aforementioned processes become increasingly dominant toward larger, more complex systems. More complex models however also require more detailed and therefore more location specific inputs, such as spatiotemporal patterns in suspended sediment concentration and seasonal effects on seagrass biomass and light requirement, which are not always available. The current model also does not account for the effects of interannual and de Smit et al.
Seagrass resilience to changing winds seasonal variability in turbidity pressure on Z. marina persistence, which is affected both by seasonal timing of winds and changes in Z. marina biomass (Hansen and Reidenbach 2013). This would however require full parameterization of disturbance-recovery timescales, and of seasonal variation in Z. marina biomass and light requirement. Hence, current model results aim to provide a general sensitivity analysis of Z. marina meadows to climate change-induced changes in wind patterns, rather than a detailed quantification of the persistence of a specific Z. marina meadow under these conditions. The results from the climate change scenarios indicate that sheltered seagrass meadows such as Bokevik bay are more sensitive to changes in average conditions, i.e., median winds, than changes in storminess in terms of SSL-feedback mechanisms. Z. marina is able to maintain its biomass up to a u cr exceedance percentage between 25% and 50%, and mortality only occurs between 50% and 75% u cr exceedance time (Biber et al. 2009). The mortality threshold is not reached in any of the climate change scenarios. In addition, the increase in P (u > u cr ) due to increasing mean wind velocity varies little with Z. marina blade area. Thus, the Z. marina meadow at Bokevik bay will not collapse or become bistable due to potential future changes in the SSL-feedback. Hence, assessing the sensitivity of SSL feedback dynamics to changing wind conditions may yield important insights in future seagrass stability even in wind sheltered systems.
Fully exposed bays without topographic sheltering showed to be more sensitive to changing wind conditions. While the effect of increasing storminess in the fully exposed Z. marina meadow with similar initial P(u > u cr ) as Bokevik bay (C 1.9 ) was still limited, increasing mean wind velocity caused an exponential increase in P(u > u cr ) toward lower blade area and bare sediment. Hence, increasing mean wind velocity leads to an increasing difference in P(u > u cr ) between bare and vegetated sediments, and may thus lead to SSL feedback induced bistability. However, the limited changes in turbidity pressure inside the Z. marina meadow indicate that existing meadows under low turbidity pressure are resilient to changing wind conditions. When Z. marina is under an initially high P (u > u cr ) (C 6 ), the increase in P(u > u cr ) was consistent across the blade area gradient, indicating that seagrass meadows under higher turbidity pressure are also highly sensitive to changes in wind conditions. While the effect of increasing storminess on Z. marina stability is limited in terms of SSL feedback mechanisms, it may still have important consequences due to its effect on other processes. For example, increasing frequency of storm events causing loss of aboveground biomass may lead to collapse if recovery time is insufficient ). Increasing storm magnitude may lead to the exceedance of other sediment-related seagrass survival thresholds such as burial (Cabaço et al. 2008), uprooting (Infantes et al. 2011), or edge scouring and cliff erosion due to increasing maximum wave energy (Chen et al. 2012;Twomey et al. 2021). As these thresholds, in contrast to the SSL feedback, are instantaneous, large individual storms can cause major mortality events, which have important implications on long-term seagrass survival (Gera et al. 2014;Oprandi et al. 2020).

Outlook and management implications
In this study, we showed that the amount of bed protection by flexible submerged vegetation across a range of densities can be predicted with a combination of measurements of sediment erodibility with removed aboveground biomass and species-specific blade area-u cr relations. Such relations can be used to provide insight in the wave climates at which critical transitions occur in the Seagrass-Sediment-Light (SSL) feedback, which induces alternative stable states in seagrass meadows (van der Heide et al. 2007;Carr et al. 2010;Adams et al. 2016Adams et al. , 2018. Hydrodynamics, seagrass morphology, seagrass density, and sediment characteristics are highly variable between sites and within seagrass meadows. Therefore, in situ measurements of u cr using field flumes are needed to identify these location specific hydrodynamic thresholds for seagrass meadow stability. For seagrass preservation, measured sediments resuspension thresholds in combination with wave predictions based on future changes in sea level and wind patterns can be used to assess to which extent existing seagrass meadows are vulnerable to changing wind conditions due to climate change. For seagrass restoration, measured sediment resuspension thresholds may be used to identify the most suitable locations for restoration efforts.

Data availability statement
The data and scripts used to generate the presented results are publicly available at the 4TU repository, doi:10.4121/ 14546124.