# Flexural Rigidity and Shoot Reconfiguration Determine Wake Length Behind Saltmarsh Vegetation Patches

Tim I. Marjoribanks, Dimitri Lague, R.J. Hardy, R.J. Boothroyd, J. Leroux, Cendrine Mony, Sara Puijalon

► To cite this version:

Tim I. Marjoribanks, Dimitri Lague, R.J. Hardy, R.J. Boothroyd, J. Leroux, et al., Flexural Rigidity and Shoot Reconfiguration Determine Wake Length Behind Saltmarsh Vegetation Patches. Journal of Geophysical Research: Earth Surface, 2019, 124 (8), pp.2176-2196. 10.1029/2019JF005012. hal-02308542

HAL Id: hal-02308542

https://univ-rennes.hal.science/hal-02308542

Submitted on 8 Oct 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2019JF005012

**Key Points:**

* We develop a new field-parameterized model for simulating shoot reconfiguration within vegetation patches

* Wake length is controlled by a vertical shear layer, which in turn depends on shoot flexural rigidity and flow velocity

* A simplified beam equation model allows prediction of wake length from velocity and flexural rigidity

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

T. I. Marjoribanks,
t.i.marjoribanks@lboro.ac.uk

**Citation:**

Marjoribanks, T. I., Lague, D., Hardy, R. J., Boothroyd, R. J., Leroux, J., Mony, C., & Puijalon, S. (2019). Flexural rigidity and shoot reconfiguration determine wake length behind saltmarsh vegetation patches. *Journal of Geophysical Research: Earth Surface, 124*, 2176–2196. https://doi.org/10.1029/2019JF005012

Received 17 JAN 2019

Accepted 2 AUG 2019

Accepted article online 6 AUG 2019

Published online 17 AUG 2019

©2019. The Authors.
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.

<figure><img src="image_2.png" /></figure>

# Flexural Rigidity and Shoot Reconfiguration Determine Wake Length Behind Saltmarsh Vegetation Patches

T. I. Marjoribanks¹ ID, D. Lague² ID, R. J. Hardy³ ID, R. J. Boothroyd⁴, J. Leroux²,⁵, C. Mony⁶, and S. Puijalon⁷

¹School of Architecture, Civil and Building Engineering, Loughborough University, Loughborough, UK, ²Université de Rennes, CNRS, Geosciences, UMR 6118, Rennes, France, ³Department of Geography, Durham University, Durham, UK, ⁴School of Geographical and Earth Sciences, University of Glasgow, Glasgow, UK, ⁵Now at Azur Drones, Boulogne-Billancourt, France, ⁶Université de Rennes, CNRS, ECOBIO, UMR 6553, Rennes, France, ⁷Université de Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR5023 LEHNA, Villeurbanne, France

**Abstract** Vegetation patches play an important role in controlling sediment deposition in shallow aquatic environments such as coastal saltmarshes and fluvial systems. However, predicting deposition around vegetation patches is difficult due to the complexity of patch morphology and their dynamic interaction with the flow. Here we incorporate a biomechanical model, parameterized using field data, within a 3-D computational fluid dynamics model which allows prediction of individual shoot reconfiguration within patches due to flow forcing. The model predicts velocity attenuation and bed shear stresses within the wake of the patch which agree spatially with accretion patterns measured in the field using terrestrial lidar. The model is applied to sparse patches of *Suaeda maritima*, located in saltmarshes of coastal habitats, to explore the role of (I) shoot distribution, (II) patch geometry, (III) shoot flexural rigidity, and (IV) bulk flow velocity in determining the length of the predicted wake region. We demonstrate that for *Suaeda maritima*, with intermediate rigidity, the vertical shear layer over the vegetation controls the length of the predicted wake region. Consequently, reconfiguration due to flexural rigidity strongly impacts on wake length, confounding the relationship between patch height and wake length. A simplified model for predicting wake length based on shoot reconfiguration is applied to the simulation data and shows good agreement. The results demonstrate that the observed wake characteristics can be well explained by intraspecific variability in flexural rigidity, thus demonstrating the importance of biomechanical traits in determining flow-vegetation-sediment interactions.

## 1. Introduction

Vegetation has a major impact on flow, sediment dynamics, and ecological processes within shallow aquatic environments and can control the evolution of fluvial, deltaic, and coastal sedimentary systems including coastal saltmarshes (Neumeier & Amos, 2006; Temmerman et al., 2005). Understanding the role of vegetation in saltmarsh development and evolution is crucial for managing natural sea defenses and mitigating against future sea level rise (Möller, 2006; Shepard et al., 2011; Temmerman et al., 2013). Saltmarsh vegetation communities develop through a process of colonization by pioneer species, which alter flow and sedimentation for succession to occur. Such sedimentation occurs behind isolated vegetation patches and controls the process of succession, enabling the establishment of species more sensitive to higher inundation frequency and shear stresses through facilitation interaction processes (Castellanos et al., 1994; Langlois et al., 2003; Tempest et al., 2015).

Previous research has shown that sediment deposition behind idealized vegetation patches may be explained using the hydraulic principle of flow separation and reattachment (Zong & Nepf, 2012; see Figure 1). Vegetation represents a porous blockage, which extracts energy from the flow as it crosses the patch, causing a wake to form behind, characterized by low bleed velocities and reduced turbulent kinetic energy (Nepf, 1999; Nicolle & Eames, 2011). Here, enhanced deposition may occur (e.g., Figure 2), depending upon local sediment and flow conditions (Chen et al., 2012; Shi et al., 2016; Zong & Nepf, 2012). Where the flow reattaches, turbulent kinetic energy is higher, and vortex shedding may occur, which limits sediment deposition (Chen et al., 2012; Zong & Nepf, 2012). Previous studies have defined the spatial extent of the wake sedimentation region using threshold values of both velocity and turbulent kinetic energy (Chen et al., 2012; de Lima
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic diagram of (a) horizontal and (b) vertical shear layer development and their control on wake length ($L_{Z&N}$, $L_V$). In each case, the wake length is controlled by patch width (w) and deflected patch height ($h_d$), respectively. $L_{Z&N}$ is the width-dependent wake length, calculated geometrically as $L_{Z&N} = \frac{w}{2} \frac{\bar{U}}{S_{\delta} \Delta U}$ (Zong & Nepf, 2012), and $L_V$ is the equivalent height-dependent wake length calculated geometrically as</figcaption></figure>

$$L_V = x_d + h_d \frac{\bar{U}}{S_{\delta} \Delta U} \quad (\text{see section 2.2.3}). \text{ Here, } S_{\delta} \text{ is the shear layer growth rate. Representative velocity profiles are indicated by the black arrows. Downstream shear layer growth is shown by dotted lines.}$$

et al., 2015; Yang et al., 2016). To date, most research has focused on the role of the horizontal shear layer and wake structure, formed by the horizontal flow separation around the patch, in determining sediment deposition extent (Figure 1a). In these cases, for given patch porosity and flow conditions, patch width determines wake length (Zong & Nepf, 2012). However, for submerged vegetation, the vertical shear layer formed at the top of the patch (Figure 1b) may be of equal or greater importance, depending on patch geometry (Hu et al., 2018; Liu et al., 2018). If wake length is controlled by the vertical shear layer instead of the horizontal shear layer, then shoot height, rather than patch width, will determine the wake length. Moreover, if the wake is controlled by shoot height, then shoot reconfiguration, the adjustment of the shoot into a more streamlined position to reduce drag (Luhar & Nepf, 2011; Nepf & Vivoni, 2000), will also affect the wake length. Therefore, we hypothesize that flexural rigidity is an important parameter in defining wake length for flexible, submerged vegetation patches.

Testing this hypothesis requires investigating the impact of flow and patch characteristics on wake length in an environment where such parameters may be quantified and controlled. Computational fluid dynamics (CFD) modeling provides an effective tool for this as it allows control of flow, patch geometry, and biomechanical properties which may be difficult to accurately constrain within field or flume environments. This permits testing over a wide parameter space within a controlled setup. However, currently such models either allow for accurate prediction of the motion of single shoots (e.g., Dijkstra & Uittenbogaard, 2010; Marjoribanks et al., 2014) or rely on a broader-scale representation of vegetation using a field-parameterized rigid cylinder drag force approach (Baptist et al., 2007; Temmerman et al., 2005; Tempest et al., 2015; Wu et al., 2017). Hence, there is a need for models which are field parameterized but also allow for prediction of flow-vegetation interactions at the shoot scale. Importantly, existing models of flow and deposition processes in the wake

of vegetation patches cannot currently account for the impacts of shoot reconfiguration which, through altering the canopy height and shear layer characteristics, may affect patch wake length.

The aim of this paper is to test the hypothesis that flexural rigidity is important in determining wake length behind flexible, submerged saltmarsh vegetation (*Suaeda maritima*) patches. This is achieved through (i) developing a field-parameterized combined CFD-biomechanical model that can represent shoot reconfiguration; (ii) using it to investigate the relative importance of flexural rigidity in determining wake length, compared to the impact of patch geometry (patch width/height and shoot distribution) and flow conditions (velocity); and (iii) evaluating the performance of a simplified analytical model for vertical shear layers, adapted from Zong and Nepf's (2012) model for horizontal shear layers (Figure 1a) and based on flexural rigidity and bulk flow characteristics, in predicting wake lengths observed in the field data and CFD simulations.

## 2. Methods

### 2.1. Field Data Collection for Model Parameterization and Validation

Our model plant species is *Suaeda maritima*, an annual halophyte plant common in the pioneer saltmarshes of the Mont Saint-Michel Bay (France, 48.631°N, -1.504°W; Figures 2a and 2c). *Suaeda maritima* grows from June to October with a peak height reached around September to October (Tessier et al., 2000). Our study area is the seaward part of a tidal meander where sparse patches of *Suaeda maritima* growing on the inner bar are directly exposed to flow crossing the inner bar during flooding. Wake deposits of several centimeter thickness develop downstream of each patch (Figures 2b and 2c) with respect to the flood direction during the spring tides of September and October when the vegetation is at its peak development. The
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Study site and Terrestrial Laser Scanner data. (a) aerial photograph of the site in September 2010 (Institut Geographique National). (b) Map of vertical topographic change between 11 and 4 October calculated on ground points with the M3C2 algorithm (Lague et al., 2013). Green points are automatically classified vegetation patches using CANUPO (Brodu & Lague, 2012). The coordinate system indicated on the picture is used to compute the vegetation patch characteristics and wake deposit length. The location of the reference ADV is indicated (ADV<sub>ref</sub>). (c) Pictures of the site illustrating the development of wake deposits during the spring tides. (d) Detailed view of the red rectangle in b. (e) Cross-section A-A' (width of 20 mm) showing accretion patterns. The vegetation blocks laser scans immediately in front of the patch. ADV = acoustic doppler velocimeters.</figcaption></figure>

wake deposits are composed of very fine sandy to muddy sediment, with a median grain size ranging from 30 to 100 µm. Owing to the meander configuration, ebb flow during the spring tides is concentrated in the outer bend and tends to circumvent the vegetated inner bar (Leroux, 2013). Consequently, flow velocity and suspended sediment concentration (SSC) remain much lower in the studied area during ebbing compared to flooding and no wake deposits develop in the opposite direction of the patch (see section 2.1.2 and Figures 2 and 3). Hence, apart from a uniform component of accretion related to ebbing that we account
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Suspended sediment concentration and (b) velocity and depth time series over a single tide cycle for locations within a patch wake and at a reference location within the freestream flow (Figure 2b and supporting information Figures S2 and S3). Dashed line shows flow depth above $ADV_{ref}$. SSC = suspended sediment concentration; ADV = acoustic doppler velocimeters.</figcaption></figure>

for in data processing, the development of wake deposits occurs solely
during flooding, under unidirectional flow.

Here, we used a combination of high-resolution Terrestrial Laser Scanner (TLS or lidar) surveys, in situ measurement of flow velocity, SSC, and plant sampling to both parameterize and validate the numerical model.

2.1.1. 3-D TLS Data Analysis

Terrestrial laser scanning provides an ideal method for the large-scale,
high-resolution survey of temporally variable sediment and vegetation
characteristics that are required to parameterize and validate a patch-
scale numerical model. Daily TLS surveys took place during spring tides
in October 2010 using a Leica Scanstation 2 placed at three different scan-
ning positions. The setup and data processing were as used by Lague et al.
(2013) using four targets anchored on consolidated marsh for local geore-
ferencing between daily surveys. Scanning time available between tides
and difficulties in scanning in the unconsolidated muddy channel depos-
its meant that scanning positions did not fully cover the area resulting in
partial shadow effects (i.e., no data) in the area of densest vegetation
patches. The point density varied spatially between 1.0 × 10³ and 1.0 ×
10⁵ points/m² (Figure 2).

We used the first survey on 4 October before the first spring tides to extract
vegetation characteristics from the point cloud using various algorithms
implemented in the open source 3-D point cloud processing software

Cloudcompare V2.5 (General Public Licence, 2014). We used the CANUPO classification algorithm (Brodu & Lague, 2012) to automatically separate vegetation patches from the ground. The CANUPO algorithm is a supervised semantic 3-D classification method operating directly on irregular point clouds. It computes the 3-D geometry of the point cloud at various scales using a principal component analysis of the points coordinates. Once trained with manually selected samples, it results in a classification for each point associated with a metric of classification confidence. On the studied vegetation patches, we used a range of scales from 3 to 6 cm with a 0.5-cm interval. This resulted in classification accuracy better than 99% (Brodu & Lague, 2012). The classified vegetation point cloud was subsequently segmented into individual patches using a connected component algorithm that groups subsets of the 3-D point cloud which have a nearest point distance smaller than a given threshold (5 cm for this data set). This resulted in 216 individual vegetation patches (Figure 2b). We assumed a uniform flow direction for the study site and considered a fixed coordinate system in which the y axis was parallel to flow direction computed as the mean direction of the tail deposits (Figure 2b). Patch width *w*, length *l*, and height *h* were measured along the *x*, *y*, and *z* directions to ensure a consistent approach for all vegetation patches with respect to the mean flow direction. The size characteristics *w*, *l*, and *h* correspond to the range between the 2nd and 98th percentiles of the coordinates (*x*, *y*, *z*) of all the points of a given patch. Using percentiles reduces the dependency on the occasional single shoot sticking out of the vegetation patch and better describes the geometry of the densest part. The observed patches of *Suaeda maritima* tended to develop mostly perpendicular to the flow direction (Figure 2c). Consequently, the width was generally larger than the length.

Wake sediment characteristics were extracted from the topographic change between 4 and 11 October of the ground points (supporting information Figure S1). The topographic change was measured vertically using the M3C2 algorithm dedicated to accurate change detection on 3-D point clouds (Lague et al., 2013). M3C2 operates directly on the two 3-D point clouds by averaging the point cloud position within a disk of specified diameter (chosen here as 20 cm) resulting in a map of topographic change (Figure 2b). By looking at topographic change we can define the wake as the zone where flow velocity and turbulent kinetic energy were low enough for sediment to settle. The map of topographic change shows that the ambient accretion rate away from wake deposits slightly varies along the flood direction such that a uniform threshold of accretion cannot simply be used (Figures 2d and 2e; Leroux, 2013). To account for the spatial variation of the ambient accretion rate, we first conducted a manual, coarse removal of the wake deposits highlighted by the M3C2 measurement on the 11 October data (Figure 2b). We used the resulting point cloud to
interpolate the 11 October ground data in the location of wake deposits. We used M3C2 with a projection scale of 20 cm to measure the distance between the 11 October ground data and the interpolated point cloud. A distance threshold of 3 mm was chosen as it best distinguished the wake deposits behind vegetation patches. Due to the ground-based viewpoint of the TLS data, the vegetation patches generated significant shadow in the densest patch area and not all wake deposits could be properly detected. We thus manually selected wake deposits and associated vegetation patches in the sparsest area ensuring that the wake deposits could be clearly separated in the case of two nearby vegetation patches. This choice is consistent with the setup of the numerical simulation which explores a single patch. This restricted the analysis to 45 vegetation patches and their associated wake deposits (supporting information Figure S1). Mean accretion on the chosen wake deposits was 40 mm over the survey period with local peak accretion up to ~80 mm (Figures 2d and 2e). The length of each patch was measured along the y axis from the center of the corresponding vegetation patch. The mean accretion rate outside the wake deposits was 20 mm.

### 2.1.2. Flow and Sediment Measurement

The field data were originally collected as part of a separate project, and therefore velocity data were not collected with detailed numerical modeling in mind. Consequently, flow and sediment measurements within the study area were limited. Two 6-MHz Nortek Vector acoustic doppler velocimeters (ADV) measuring the three components of the velocity field were installed between 6 and 10 October and set up to record 1-min burst measurements at 64 Hz every 10 min. One ADV was installed immediately behind a very large patch of vegetation (ADV<sub>wake</sub>), and the other directly exposed to the incoming flood and away from the influence of any vegetation wake (ADV<sub>ref</sub>; supporting information Figures S1 and S2). The ADVs were fixed on rigid frames and set up to measure the velocity at 0.06 m above ground. After the survey period, this distance was 0.05 m for ADV<sub>ref</sub> and 0.005 m for ADV<sub>wake</sub> illustrating the large difference in accretion between the two locations. ADV raw data were filtered to remove low signal to noise ratio data. SSC was estimated using the ADVs backscatter signal averaged over each 1-min burst measurement (Salehi & Strom, 2011) and previously calibrated against water sampling during one tide event (Leroux, 2013). Figure 3 shows a typical record during a tidal cycle. Peak velocity (~0.6 m/s) and SSC (~2.8 g/L) away from the wake of vegetation occurs during flooding. At the same time, the SSC in the wake zone is slightly reduced (~1 g/L), and the flow velocity strongly damped (~0.1–0.15 m/s). The combination of reduced velocity and high SSC is conducive of high accretion rates in the wake zone. During ebbing, the two ADVs exhibit similar characteristics with reduced peak flow velocity (~0.2 m/s) and a much lower SSC (~0.2–0.4 g/L). Leroux (2013) showed that the peak flood velocity away from the vegetation wake increases linearly with the tide amplitude (supporting information Figure S3), while the peak flood velocity in the vegetation wake is nearly constant, around 0.15 m/s.

### 2.1.3. Stem Geometrical and Mechanical Properties

As the biomechanical analysis of the plants was not initially planned when the TLS surveys started in 2010, it occurred 2 years later in July 2012 on vegetation patches located on the same area. The density of vegetation patches was higher, but the samples were chosen in the more exposed and sparse area to be as close as possible to the configuration of October 2010. The flexural rigidity of 196 stems sampled in two 3.5 × 3.5-m² squares was measured with a universal testing machine (Instron 5942 Canton, MA, USA). Basal stem samples (50 mm long) were tested as cantilever beams (one fixed end bending test; Hamann & Puijalon, 2013). Each sample was clamped horizontally at its basal end while a force was applied at its midpoint by lowering a probe at a rate of 10 mm/min. The following biomechanical traits were calculated:

1. The bending Young's modulus, *E* (Pa), quantifies the material stiffness and is calculated as the slope of the stress-strain curve in the elastic deformation region;

2. The second moment of area, *I* (m⁴), accounts for the effect of the cross-sectional geometry of a structure on its bending stress. As stem cross section was circular, *I* was calculated as *I* = (π*r*⁴)/4, where *r* is the radius of stem cross section (Niklas, 1992);

3. The flexural rigidity, *EI* (Nm²), quantifies the stiffness of the stem fragment and was calculated by multiplying *E* and *I*.

A second sampling took place in September 2013 to evaluate the shoot density and length in a vegetation state more similar to that when measurements were made in October 2010 and to constrain the equivalent porosity of the shoots for the numerical modeling setup.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Relationship between (a) patch length and width and (b) patch height and width for *Suaeda maritima* patches at Mont Saint Michel. The black lines indicate regression line fits for the equations given.</figcaption></figure>

## 2.2. CFD Modeling Methods

### 2.2.1. Vegetation Model Conceptualization and Parameterization From Field Data

Each vegetation patch was represented by a given planform area, scaled according to the field data. The patch characteristics measured from the TLS data (Figure 4) show a linear relationship between patch width and patch length, such that $l = 0.6w$ ($r^2 = 0.87$), while the relationship between the patch height and width exhibits a weaker power law trend, $h = 0.25w^{0.6}$ ($r^2 = 0.59$). Although there is some scatter, the data capture the mean geometrical characteristics of the patches. From the TLS data, these characteristics appear independent of the density of vegetation patches.

For each simulation, the given patch area ($l \times w$, see Figure 5a) was populated with individual shoots matching the mean field-observed shoot density of 505 ± 202 shoots/m² measured over 15 patches of width greater than 0.1 m. The location of each shoot within the patch area was randomly selected from a uniform distribution in both the streamwise and spanwise directions, with the additional criteria that shoots could not overlap and must be separated in each direction by 5 mm (one numerical grid cell) of unoccupied space. Shoot properties were assigned using the characteristics observed in the field (Figure 5b). For shoot height, an average spanwise slope of 0.26 in shoot height (Δ Shoot height) was observed away from the patch center and therefore this was applied within the model such that $h = 0.25w^{0.6}$ was the maximum shoot height at the patch center. The stem flexural rigidity, calculated using samples from the field site, was $EI=1.56 \times 10^{-4} \pm 1.5 \times 10^{-4} \text{ Nm}^2$. However, the distribution was not Gaussian and spanned approximately 2 orders of magnitude from $1.0 \times 10^{-5}$ to $1.0 \times 10^{-3} \text{ Nm}^2$.

Each vegetation shoot was represented within the model as a dynamic porous blockage which reconfigured to the mean flow at each model iteration. Within each shoot, the individual main stem and leaf features were

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Incorporation of vegetation patches within the numerical model: (a) Average patch dimensions obtained from the laser scan field data. (b) Further patch parameters were obtained from both laser scan data and photographs and used to generate random, statistically similar model patches that were incorporated into the CFD model (section 2.2.2). (c) The numerical model domain, which was scaled by patch width ($w$) and contained a single patch of shoots. The example shoot shows the arrangement of shoot centered (1), intermediate (2), and shoot edge (3) porosity zones within each shoot. Diagram not drawn to scale. CFD = computational fluid dynamics.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Image analysis process applied to sample shoots. (a) Manual identification of shoot main stems (dashed lines). (b) Binary image created using RGB masking. (c) Solid volume fraction map. (d) Extracted shoots with curvature removed. (e) Resampled solid volume fraction at 5-mm resolution. (f) Mean solid volume fraction at each location within the shoot. (g) Standard deviation of solid volume fraction at each location within the shoot.</figcaption></figure>

predominantly subgrid scale and therefore their effects on the flow were modeled using a combined dynamic mass flux scaling algorithm (Hardy et al., 2005; Marjoribanks et al., 2014) and an additional drag force term within the Navier-Stokes equations which represented the mass and momentum effects of the vegetation respectively. Each vegetation shoot was conceptualized as a vertical stack of porous arrays, 5 × 5 numerical cells wide (0.025 × 0.025 m, 5Δx × 5Δx; see Figure 5c). This width was chosen based on the average width measured from plants sampled in the field. The center of each array corresponded to a single node on a beam which was used to calculate shoot reconfiguration and was fixed to the bed. Static reconfiguration of each shoot to the time-averaged flow was calculated by solving the steady Euler-Bernoulli beam equation (Marjoribanks et al., 2014), applying the total force acting over each porous array at each single vertical node. A relaxed iterative procedure, whereby shoot posture was updated at every flow calculation iteration, was used to ensure that the shoot positions converged simultaneously to the flow. As shoot position was directly coupled with the flow, the flow convergence criterion was sufficient to imply shoot position convergence. Shoot collisions were not treated within the model.

In order to parameterize the shoot porosity ($\phi$) for the mass flux scaling algorithm, images of three characteristic shoots, taken in the field, were analyzed as follows (Figure 6a). Within the image, the main stem of
each shoot was manually identified and interpolated into a continuous profile. The image was then converted into a binary mask using RGB thresholding (Figure 6b; Boothroyd et al., 2017). This binary mask was converted into a shoot solid volume fraction ($1 - \phi$) map (Figure 6c) for each pixel, by sampling over a 0.005-m ($\Delta x$) moving window. Using this map, a vertical profile for each shoot was created (Figure 6d), by extracting a 0.025-m horizontal window centered on the interpolated main stem at each pixel up-shoot. The results of this process (Figure 6d) visually demonstrate that using a window width of 0.025 m captures the majority of each shoot with minimal overlap. Nevertheless, there is some inherent error in the process due to the close proximity of the shoots in the original photograph. The vertical profile was then resampled horizontally into five cells, each 0.005 m wide ($\Delta x$; Figure 6e). Using these profiles, the mean (Figure 6f) and standard deviation (Figure 6g) of the shoot solid volume fraction were calculated over all three shoots for the three different relative cell positions shown in Figure 5c, representing grid cells at the shoot center (1), intermediate grid cells (2), and grid cells at the shoot boundary (3). Each porous array was randomly populated according to its position relative to the shoot center assuming a normal distribution of porosity with position-dependent mean and standard deviation values. The porosity ($\phi$) was applied at grid cell faces as outlined in Marjoribanks, Hardy, Lane, and Tancock (2017).

The drag force acting within each cell was calculated based upon the simplified assumption that all plant components were cylindrical and therefore for each the drag coefficient ($C_D$) was equal to 1 (Panton, 1996). Suaeda maritima typically exhibit a range of leaf forms from semiterete to terete (cylindrical; Tison & de Foucault, 2014), and Figure 6a shows that many of the stems and leaves were approximately cylindrical in shape. Drag coefficient will vary between stems/leaves and with orientation, and so choice of drag coefficient may impact upon the accuracy of results. However, through explicitly representing the large-scale impacts of shoot reconfiguration and morphology, we have reduced the impact of drag coefficient to the subgrid leaf and stem scale, compared to other methods where all these effects are represented through a single drag coefficient.

The drag force within every grid cell was calculated by summing up the drag due to each cylindrical stem and leaf within the cell, assuming stems and leaves were perpendicular to the flow:

$$F_D = n0.5\rho C_D h d U^2 \quad (1)$$

where n is the number of stems and leaves, $\rho$ (kg/m³) is the fluid density, $C_D$ is the drag coefficient, h (m) is the stem height, d (m) is the stem diameter, and U (m/s) is the local grid cell velocity.

In order to relate the drag force to porosity, four further simplifying assumptions were made about the stems and leaves. First, it was assumed that stems and leaves had the same diameter (d). Second, it was assumed that they extended across either the full width ($h = \Delta y$) or height ($h = \Delta z$) of a grid cell, depending on orientation. Third, it was assumed that stems and leaves within a grid cell did not overlap, and fourth, it was assumed that they did not induce any within-cell sheltering effects. This allowed calculation of the number of stems per grid cell volume from the cell face porosity ($\phi$) as follows:

$$n = \frac{(1 - \phi) \Delta y}{d} \quad (2)$$

Here $\Delta y$ (m) could be used interchangeably with $\Delta z$ ($\Delta y = \Delta z$) depending on stem/leaf orientation. Therefore, combining equations (1) and (2), the drag force within each grid cell could be expressed as follows:

$$F_D = 0.5\rho C_D (1 - \phi) (\Delta y)^2 U^2 \quad (3)$$

This additional drag term was incorporated within the Navier Stokes equations using a linearized source term approach as outlined by Marjoribanks, Hardy, Lane, and Tancock (2017).

### 2.2.2. Flow Model

Flow was modeled using the Reynolds-averaged Navier-Stokes (RANS) equations with a Re-Normalisation Group theory derived $k - \epsilon$ turbulence closure model to produce a time-averaged, steady state solution. The RANS equations were solved using a staggered grid and a hybrid-upwind differencing scheme. The mass and momentum equations were coupled using the SIMPLEST algorithm (Spalding, 1980), and the solution was
solved iteratively until errors were reduced to 0.1% of the inlet flux. The numerical domain was scaled according to patch width to ensure comparability between scenarios, such that the domain spanned 24w in the downstream direction, 5w in the spanwise direction, and 0.75w^0.6 (3h) in the vertical direction to minimize and standardize blockage ratio and boundary effects. Thus, the blockage ratio of the vegetation patch was 0.2, though this does not consider the porosity of the patches themselves. The vegetation patch was centered at x = 2w, y = 2.5w, where in contrast to the TLS processing, the conventional frame of reference with x as the downstream coordinate and y as the cross-stream coordinate was used. For the section of the domain containing the patch, a regular grid of resolution Δx = Δy = Δz = 0.005 m was used. Downstream of the wake region (x > 4w), the x direction grid spacing was increased gradually for computational efficiency. Thus, the total number of grid cells for each simulation was between 0.46 and 38.6 × 10^6 and was 7.45 × 10^6 for the standard (w = 0.3 m) case. The computational time for each simulation scaled with the number of grid cells and was 30 hr for the w = 0.3-m case. A standard rigid lid treatment was applied at the free surface, while the bed was represented using a no-slip boundary condition and logarithmic wall function (dimensionless wall distance, y_+ > 34 for all cases). The sides of the domain were treated as frictionless boundaries. A constant inlet velocity was applied as only single-point velocity field data were available, and the outlet was represented as a fixed pressure boundary.

### 2.2.3. Wake Deposition Length Measurement

The field data presented here consist predominantly of high-resolution sediment deposition records, and therefore in order to compare these with the numerical flow data, we calculated wake lengths which we suggest correspond to potential regions of enhanced deposition. In this study we used three different measures of wake length. Primarily, wake deposition was related to areas of the bed where the bed shear stress did not exceed the critical bed shear stress, τ_c (Pa). The bed shear stress, τ_0 (Pa), was calculated using the turbulent kinetic energy per unit mass k (m^2/s^2) and may be written as

$$ \tau_0 = C_1 \rho k \quad (4) $$

The turbulent kinetic energy is a function of the fluctuating velocity components (u', v', w') and is time averaged such that $k = \frac{1}{2}(\overline{u'^2} + \overline{v'^2} + \overline{w'^2})$. It is calculated directly by the Re-Normalisation Group $k - \epsilon$ turbulence closure model by solving a transport equation, and the near-bed value was evaluated in the first grid cell above the bed (z = 0.0025 m). Recent studies have shown that near-bed turbulent kinetic energy is a good predictor of sediment transport within vegetated channels (Tinoco & Coco, 2018; Yang & Nepf, 2018). Previous field experiments within oceanic environments have found C_1 = 0.19 - 0.2 (Pope et al., 2006; Soulsby, 1981), whereas a theoretical derivation, assuming approximate local equilibrium between turbulence production and dissipation, gives a value of C_1 = 0.3 (Baranya et al., 2012; Rodi, 1980). Here we used C_1 = 0.3 in line with the CFD model wall parameterization. The critical shear stress was τ_c = 0.131 N/m^2 and was chosen to represent the threshold for motion for the largest grain size fraction (0.1 mm) applying the adapted Shields (1936) method of Soulsby (1997). Wake length, L_TKE, was determined by the furthest point downstream where the shear stress criterion was met. Sediment cohesion was not considered here as the model was used to predict reentrainment of mobile grains rather than erosion.

In addition to the shear stress method, wake length was also calculated using the theoretical basis illustrated in Figure 1a. Using this method, for horizontal shear layers, wake length may be calculated theoretically as the convergence point of the two shear layers formed at the edge of the patch (Chen et al., 2012; Zong & Nepf, 2012).

$$ L_{Z\&N} = \frac{w/2}{S_{\delta}} \frac{\bar{U}}{\Delta U} \quad (5) $$

Here w (m) is the patch width, U (m/s) is the arithmetic mean velocity of the mixing layer, ΔU (m/s) is the velocity difference across the mixing layer, and S_δ is the rate of growth of the mixing layer half-width, as measured from the starting location of the shear layer (i.e., any movement of the center of the shear layer is represented within the growth rate). However, considering the vertical shear layer (Figure 1b), equation (5) can be altered to calculate the distance downstream at which the vertical shear layer reaches the bed.
$$L_V = x_d + h_d \frac{\bar{U}}{S_\delta \Delta U} \qquad (6)$$

The second term in equation (6) is directly analogous to the Zong and Nepf (2012) model, with patch half-width replaced by the deflected shoot height, $h_d$ (m), while the first term $x_d$ (m) accounts for the fact that the downstream location of the shear layer will alter for flexible vegetation (see Figure 1b). This is similar to the term added by Hu et al. (2018) to account for displacement of the wake limit due to leaf length. For any given patch, both $L_{Z&N}$ and $L_V$ can be calculated and comparing the two predicted wake lengths enables classification of each patch wake length as either width dependent ($L_V > L_{Z&N}$) or height dependent ($L_V < L_{Z&N}$). In other words, whichever shear layer converges furthest upstream will limit the length of the wake. In order to compute $L_V$ and $L_{Z&N}$, values of $S_\delta$ were calculated for both the vertical and horizontal shear layers, along with the case-specific $\bar{U}$ and $\Delta U$ values. Shoot height varied within each simulated patch, so the maximum value within the patch was used as this has been shown to represent the effective canopy height (Hamed et al., 2017). To calculate $S_\delta$, the shear layer growth rate, the width of the shear layers ($\delta_x, \delta_z$) were calculated at 0.05-m intervals downstream of the patch. The edge of the shear layer was measured quantitatively as the point $z$ (or $y$ for lateral shear layers), at which $U(z) = U_1+0.1\Delta U$, where $U_1$ is the mean velocity of the slower of the two coflowing regions forming the mixing layer (Pope, 2000; Zong & Nepf, 2012). This is a standard method for calculating shear layer width but is not directly related to a specific turbulent kinetic energy threshold. Therefore, it is possible that the spread of the critical threshold in turbulent kinetic energy, which in turn affects the critical shear stress, may occur at a smaller or larger rate.

As the simulations were steady state and run at a constant velocity, while the field-site experienced tidal variations, an important consideration was what model velocity should be applied to most effectively reproduce the observed deposition patterns. The study area is a depositional environment, which experiences predominantly unidirectional but varying flow conditions and which exhibits net deposition over individual tide cycles (Leroux, 2013). Therefore, under many velocity conditions, wake sedimentation may not be distinguishable from deposition on the bare unvegetated bed. This phenomenon was studied by Shi et al. (2016), who found both sediment and flow controls on preferential wake deposition. Determining a threshold condition for preferential wake deposition is key to comparing predicted wake deposition extents with field data. SSC field data for vegetated and bare areas provided some detail on the threshold for preferential wake deposition. Figure 3 shows an example of the variation in SSC over an individual tide cycle. The largest difference in SSC between the patch and freestream flow, which we suggest relates to the period of most preferential deposition, occurs at the highest velocities (0.4–0.6 m/s). Below 0.4 m/s, at the end of the tide cycle, the SSC converges. However, the velocity was measured close to the bed ($z = 0.06$ m) and therefore likely underestimates the mean velocity. Therefore, the primary velocity for investigation was chosen as 0.7 m/s. However, other velocities (0.2–0.6 m/s) were also used for comparison. A single velocity condition, rather than an entire tide cycle, was selected to better highlight, and isolate, the role of velocity and other variables in controlling wake length.

### 2.2.4. Summary of CFD Simulations

To achieve the aim of this study, simulations were conducted with varying patch (width, shoot configuration, and flexural rigidity) and flow (velocity) conditions. When varying each of the different controlling variables, all other parameters were kept constant. The baseline comparison values used were a patch width of 0.3 m, a shoot flexural rigidity of $1.5 \times 10^{-4} \text{Nm}^2$, and a velocity of 0.7 m/s. The baseline width and flexural rigidity values correspond to the mean values for the field data, while the velocity was chosen as it corresponds approximately to the maximum preferential deposition condition as explained in section 2.2.3. Since the vegetation model randomly assigns shoot positions within each patch, when comparing patches with different widths and therefore different shoot distributions, it is important to understand the influence of subpatch shoot distribution on wake flow patterns which may confound any relationship with patch characteristics. Therefore, 10 different replicate simulations were conducted using identical baseline biomechanical, patch geometry, and velocity conditions with different random shoot distributions to quantify the impact of shoot distribution on wake length. For all other simulations with a patch width of 0.3 m, an identical shoot distribution was used. In total, results from 27 different simulations are presented here, split into four different comparison groups, reflecting investigation of different controlling variables: (I) shoot
distribution effects, (II) width effects, (III) flexural rigidity effects, and (IV) velocity effects (supporting information Table S1).

## 2.3. Simplified Analytical Model

For patches with horizontal shear layers, Zong and Nepf's (2012) method provides a simple model for predicting wake length. Assuming instead that wake length is controlled by the vertical shear layer and consequently is a function of shoot height, then equation (6) can be used as an equivalent simplified model for predicting wake length. The key difference between this model and that of Zong and Nepf is that, in contrast to patch width which will remain constant for vertically aligned shoots, in the case of the vertical shear layer, both $x_d$ and $h_d$ will vary with hydrodynamic forcing and flexural rigidity. By simplifying patches to a single shoot with identical length and biomechanical properties to that used in the CFD simulations, a steady state Euler-Bernoulli beam equation was used to calculate $h_d$ and $x_d$ and consequently predict wake length using equation (6).

For this simplified model, the inlet velocity was constant with depth and set equal to the CFD inlet condition ($U = 0.7$ m/s). Similarly, the model used a constant expression for drag force along the shoot, based upon the mean inlet velocity.

$$F_D = 0.5\rho C_D A U^2 \quad (7)$$

The drag coefficient was taken as $C_D = 1$, while the area of each section of the shoot, $A$ (m²), was adjusted with shoot inclination to represent the projected shoot frontal area. As with the CFD model, there are limitations with assigning a drag coefficient equal to 1; nevertheless, for the simplified model, while the drag coefficient may affect the exact wake lengths calculated, it should not substantially affect the trends produced by the model for varying rigidity. The spatial resolution of the simplified model was set equal to that used in the CFD model ($\Delta z = 0.005$ m). The computational time for the simplified model under baseline conditions was 0.15 s.

As discussed in section 2.2.3, the representative shear layer growth rate $S_\delta$ for calculating wake length using equation (6) will vary depending on the threshold turbulent kinetic energy and consequently bed shear stress. Therefore, different values of $S_\delta$ notionally represent different bed shear stress thresholds and consequently different grain size erosion thresholds. Accordingly, sensitivity to this parameter was examined through using four different values for $S_\delta$. Values of $S_\delta$ were chosen based on the relationship between $L_{TKE}$ and $L_V$ within the CFD simulations in Group I. Average values of $\overline{U}/\Delta U$ at each rigidity were obtained a priori from the CFD simulations with values varying between 0.77 and 1.03. However, using a single averaged value of 0.98 produced a similar result. The output from the simplified model was compared to wake lengths ($L_{TKE}$) calculated from the CFD model. In addition, the model was applied to the field data using the observed flexural rigidity value. However, the field data for flexural rigidity showed considerable uncertainty. Therefore, we also used the model to calculate the expected flexural rigidity value for each patch given its observed wake length. This was done by running the model over a wide parameter range for flexural rigidity (0.01–9.5 × 10⁻⁴ Nm²) for each patch and selecting the flexural rigidity with the lowest wake length error.

# 3. Results

## 3.1. CFD Model Results

### 3.1.1. General Description of the Flow

Simulation results (Figure 7) demonstrate the typical patch layout and shoot reconfiguration across a range of flexural rigidities for the baseline geometry and velocity conditions ($w=0.3$ m, $U=0.7$ m/s). Shoots reconfigure to the time-averaged flow, due to the drag force acting on the patch, demonstrated by the pressure gradient across the patches (Figures 7a–7c). Above the patch, horizontal and vertical shear layers form, which grow downstream, as demonstrated by the turbulent kinetic energy distribution (Figure 7d). Behind the patch, a wake forms, where flow velocity (Figure 7e) is reduced. The extent of this wake is determined by the rate of growth of the shear layers.

### 3.1.2. Comparison With Field Data

The simulations were designed to replicate the ensemble characteristics of the Mont Saint-Michel vegetation, rather than individual plant patches. Nevertheless, it is important to compare the simulation results,
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a–c) Simulation data showing shoot posture and bed pressure distributions for (a) EI=1.5 × 10⁻¹ Nm², (b) EI=1.5 × 10⁻⁴ Nm², and (c) EI=1.5 × 10⁻⁶ Nm². (d) Turbulent kinetic energy slices along the patch centerline and 0.01 m above the bed for EI=1.5 × 10⁻⁴ Nm². Calculated wake length for this case ($L_V$ = 2.15w) is shown. Shoots are represented as porosity isosurfaces ($\phi$ = 0.6). (e) Velocity magnitude slices for the same case as (d).</figcaption></figure>

both velocities and wake lengths, to see whether they fit within the observed parameter space and to validate
the assumptions used in modeling. Limited flow data from behind and adjacent to an individual patch were
available for one similar but slightly lower flow condition ($U_{z=0.06m}=0.61$ m/s). These show that observed
wake velocities were comparable (blue cross in Figure 8a) to those simulated. The shape of the
longitudinal velocity profile in Figure 8a also agrees well with previous studies of flow around porous

<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) Streamwise velocity along the patch centerline at *z* = 0.06 m (black line) and standard deviation between the five cases (gray), with vegetation patch location shown in red and field data point by the blue cross. The dotted line shows the mean *L*<sub>TKE</sub> value across the 10 replicates. (b) Variation in standard deviation in wake length with increasing numbers of replicates for the three wake length metrics. (c) Distribution of wake length values for the three wake length metrics showing the range (error bars) and mean values (crosses). TKE = turbulent kinetic energy.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Comparison of simulated wake lengths against field-measured wake lengths for different patch widths. Circles represent field patch data. Red crosses are simulated data at different patch widths with constant velocity, *U* = 0.7 m/s, blue crosses represent simulated data at *w* = 0.3 m with variable velocities (m/s) as labelled. The red error bars indicate the range of values obtained for the baseline conditions (*w* = 0.3 m, *U* = 0.7 m/s) with different random shoot distributions. (b) Field-measured relation between wake length and patch width for patches with *w* < 2h. (c) Field-measured relation between wake length and shoot height for patches with *w* > 2h.</figcaption></figure>

obstacles (e.g., Chen et al., 2012; de Lima et al., 2015) showing a decrease in velocity just in front of the patch, a large velocity deficit in the wake and gradual recovery that extends far downstream.

The simulated wake lengths for baseline conditions (Figure 9a, red crosses) lie within the observed data range across a variety of patch widths, suggesting the ensemble parameterization of patch geometry, threshold velocity, and flexural rigidity is broadly consistent. However, as discussed above it is not possible to assess the accuracy of the model for individual patch wakes.

### 3.1.3. Evaluation of Shoot Distribution Effects (Group I)

The impact of the different shoot distributions on patch centerline velocities is shown by the standard deviation in Figure 8a, which is largest in the wakes, where the shoot-scale heterogeneity introduces variations in the velocity field between replicates. The maximum standard deviation occurs immediately downstream of the patch with a value of 0.12*U*₀. However, this drops to 0.045*U*₀ by two patch widths (2w) downstream of the patch and to 0.03*U*₀ by three patch widths (3w) downstream. Comparing the trend in standard deviation over different numbers of replicate simulations (Figure 8b) demonstrates that for all three wake measures, the standard deviation levels off by 10 replicates.

The three different wake length metrics (Figure 8c) show several important features. First, the patch wake length is clearly height dependent, with *L*<sub>V</sub> < *L*<sub>Z&N</sub> for all shoot configurations. Therefore, it is most relevant to compare the theoretical value for the vertical shear layer, *L*<sub>V</sub>, with *L*<sub>TKE</sub>. There is a difference between the mean values for *L*<sub>V</sub> and *L*<sub>TKE</sub> with a ratio of approximately 1.5. The two measures also show a difference in range, *L*<sub>V</sub> having a range of 0.62w compared to 0.92w for *L*<sub>TKE</sub>. The larger range of values for *L*<sub>TKE</sub>, and larger standard deviation, may be explained by the lower mean value for wake length, meaning it occurs in a region closer to the patch where flow variability is higher (Figure 8a). Comparison with two other wake length metrics commonly used in the literature can be found in the supporting information (Text S1 and Figure S4).

### 3.1.4. Geometric (Width) Controls on Wake Length (Group II)

The previous section showed that for the baseline case, the wake length was height dependent. Here we broaden the analysis to the full range of patch widths investigated (*w* = 0.1 – 0.6 m). Considering the theoretical models for shear layer growth (Figure 1), for a given flow velocity and vegetation porosity the two geometric variables that determine wake length for rigid shoots are patch width and shoot height. From Figure 1 b it can be seen that for the vertical shear layer, the maximum half-width to which the shear layer can grow, before the shear layer interacts with the bed, is equal to the rigid shoot height, *δ*<sub>MAX</sub>(*h*) = *h*. For the horizontal shear layers, the maximum distance each shear layer can grow laterally before interfering is *δ*<sub>MAX</sub>(*w*) = *w*/2. The controlling factor for wake length will therefore depend on the ratio of *w*/*h* (Hu et al.,
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Theoretical and observed shear layer width characteristics. (a) The theoretical limit on maximum shear layer width ($\delta_{MAX}$) of the vertical and horizontal shear layers for both undeflected and deflected cases, governed by patch height ($h$), deflected patch height ($h_d$), and patch width ($w$). Note that the maximum shoot heights are used to represent the patch heights. Error bars show the range of maximum deflected shoot height between different shoot distributions. (b) Simulated vertical ($\delta_z$) and horizontal ($\delta_x$) shear layer width against distance from the end of the patch ($x_0$) for an example case. Three different theoretical growth rates ($S_\delta$) are plotted for comparison.</figcaption></figure>

2018; Liu et al., 2018). If the vertical and horizontal shear layer characteristics ($\overline{U}/S_\delta\Delta U$) are identical between the horizontal and vertical shear layers, then it follows that the smallest value of $\delta_{MAX}$ will dictate the wake length. Given the geometric relationship between $w$ and $h$ within the CFD model, based on the field data (Figure 4), for all patch widths greater than 0.177 m, the theoretical models suggest that wake length is controlled by shoot height ($\delta_{MAX}(h) < \delta_{MAX}(w)$; Figure 10a).

For flexible vegetation, the deflected shoot height replaces the rigid height ($\delta_{MAX}(h_d) = h_d$). The height dependence in Figure 10a is accentuated by shoot bending which further reduces the shoot height and therefore maximum shear layer width. When considering the simulated deflected shoot heights, $\delta_{MAX}(h_d) < \delta_{MAX}(w)$ for all patch widths simulated. The replicates for multiple shoot distributions at $w = 0.3$ m suggest that this trend is independent of shoot distribution, with a range in maximum deflected shoot height of 0.01 m between different shoot distributions. However, for deflected vegetation, equation (6) shows that the location of the shear layer will be displaced downstream ($x_d$) and therefore this will counteract the reduction in shear layer width, lengthening the wake. Therefore, the value of $x_d$ and the magnitude of the shear layer growth rate may still cause the wake to be width dominated for patches where $x_d$ is large.

The preceding analysis assumes that the vertical and horizontal shear layers have a similar growth rate, which may not necessarily be true given the difference in patch boundary roughness, foliage, and shoot reconfiguration (e.g., Hu et al., 2018). However, simulated growth of shear layer width ($\delta_x, \delta_z$), quantified using the method described in section 2.2.3, shows that for an example case, both the vertical and horizontal shear layer growth rates are similar (Figure 10b). On average, after initial variation, both can be well explained by a growth rate of $S_\delta = 0.1$. Zong and Nepf (2012) also found $S_\delta \approx 0.1$, which is within the broad range of theoretical mixing layer spreading rates observed in previous shear layer studies $S_\delta = \alpha = 0.06-0.12$ (Dimotakis, 1991; Marjoribanks et al., 2014; Pope, 2000; Sukhodolova & Sukhodolov, 2012). Therefore, it is justifiable to state that for rigid shoots, patches with $w> 0.177$ m will have height-dependent wake lengths. For flexible shoots, further analysis is required to examine the impact of reconfiguration on $h_d$ and $x_d$ and therefore wake length.

### 3.1.5. Effect of Flexural Rigidity on Wake Length (Group III)

Given the dependence of wake length on shoot deflection, it is important to understand how flexural rigidity changes deflected shoot height and consequently wake length. The CFD simulation output (Figure 7) shows the clear impact of three different flexural rigidities on patch posture, with the patch becoming increasingly streamlined with lower flexural rigidity. There is also a change in pressure distribution around the patch with larger zones of high and low pressure at the front and rear of the canopy for the more rigid cases, indicating higher patch-induced drag.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (a) Effect of flexural rigidity on maximum and mean deflected shoot height, normalized by the rigid shoot height (h). (b) Effect of flexural rigidity on wake length across five different critical shear stress (N/m²) thresholds. For τ = 0.1 N/m², there was no wake formed for EI > 3.0 × 10⁻⁴ Nm². (c) Measured distribution of flexural rigidity values at the field site. (d–f) Distribution of turbulent kinetic energy at the bed and along the patch centerline for (d) EI = 1.5 × 10⁻¹ Nm², (e) EI = 1.5 × 10⁻⁴ Nm², and (f) EI = 1.5 × 10⁻⁵ Nm².</figcaption></figure>

Simulations across eight different flexural rigidities show a clear relationship between flexural rigidity and wake length (Figure 11). As the $L_{TKE}$ definition of wake length depends directly on critical shear stress, wake lengths for six different values of critical shear stress, including the field estimated value of $\tau_c = 0.13$ N/m², are plotted for comparison (Figure 11b). The trends for all critical shear stress values are similar. At low rigidities ($EI<3 \times 10^{-4}$ Nm²), wake length increases with rigidity. A maximum wake length then occurs at approximately $EI = 3 \times 10^{-4}$ Nm², above which the wake length decreases before converging to a constant value. The wake length is constant for high rigidities ($EI > 1.5 \times 10^{-2}$ Nm²), as most shoots have reached their fully rigid positions (Figure 11d). Critical shear stress has a variable effect on wake length, depending on flexural rigidity, with greater variability evident at higher rigidities. This suggests that the longitudinal gradient of bed shear stress is steeper for highly flexible canopies.

### 3.1.6. Effect of Velocity on Wake Length (Group IV)

The results with varying velocity (Figure 12) show that velocity in the wake region, downstream of the individual shoot wakes, is relatively independent of inlet velocity, which matches the results from the limited field ADV data (supporting information Figure S3). Varying inlet velocity for a constant flexural rigidity (Figure 12) has a similar impact on patch reconfiguration and deflected patch height to the case in the previous section with constant velocity and variable rigidity (Figures 11d–11f). This is because both represent an increase in Cauchy number, which is the ratio of drag to rigidity forces and has been shown to determine shoot reconfiguration (see Luhar & Nepf, 2011). However, the impact of velocity on wake length is more complicated as variations in inlet velocity affect the bleed flow and therefore shear layer characteristics ($\Delta U, \bar{U}$). Furthermore, at low velocity cases, turbulent kinetic energy levels are sufficiently reduced (Figure 12b) such that preferential wake deposition may not occur. Wake length, $L_{TKE}$, decreases with inlet velocity (supporting information Figure S5), but this trend is not discussed further here because of the additional complexities outlined above.

## 3.2. Simplified Analytical Model Results

Based on the observed average relationship between $L_V$ and $L_{TKE}$ in the CFD simulation data (Figure 8c), within the simplified analytical model, values of $S_δ ≈ 0.15 – 0.16$ were chosen such that the wake lengths
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The impact of inlet velocity, *U*, on (a) wake velocity distribution and (b) turbulent kinetic energy distribution for three different inlet velocity conditions.</figcaption></figure>

produced by the simplified model using the *L*<sub>V</sub> method should approximately equate to *L*<sub>TKE</sub> values calculated from the CFD data. Applying the model across a range of flexural rigidities and shear layer growth rates (Figure 13a) produces a graph similar in shape to Figure 11 and predictions that visually agree well with the simulated wake lengths for the field-estimated critical shear stress, τ = 0.131 N/m² (black line in Figure 13a) particularly at lower rigidities (EI<2.0 × 10<sup>-4</sup> Nm²). At higher rigidities, the simulated decrease in wake length is underrepresented by the simplified model. The results with different shear layer spreading rates, *S*<sub>δ</sub>, show that the impact of this parameter differs with flexural rigidity, with the difference in wake length between the four cases varying between 0.1* *w* at EI = 1.5 × 10<sup>-5</sup> Nm² and 0.44* *w* at EI = 3 × 10<sup>-4</sup> Nm². The simulation data plot within the range *S*<sub>δ</sub> = 0.14–0.17 but do not follow the trend of a single value. Instead, the simulation data agree best with *S*<sub>δ</sub> = 0.15 at intermediate rigidities and *S*<sub>δ</sub> = 0.17 at both high and low rigidities.

Applying the simplified model, with $S_{\delta}$ = 0.15, to each of the 45 field-surveyed patches with an identical flex-
ural rigidity of 1.5 × 10⁻⁴ Nm² produced an average error in predicted wake length of 0.22 m, which is
expected given the local variability in both flexural rigidity and velocity. However, using the analytical model

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Simplified modeling of the relationship between flexural rigidity and wake length. (a) Comparison of predicted (colored symbols) and simulated (black line) wake length with different wake growth rates (*S*<sub>δ</sub>). (b) Schematic diagram illustrating the relative importance of variation in *x*<sub>d</sub> and *h*<sub>d</sub> in determining the change in wake length in upright and prone positions. (c) Flexural rigidities calculated from the observed patch wake lengths plotted alongside the observed distribution of flexural rigidities from the stems sampled in the field. The red dashed line shows the mean value from the field data.</figcaption></figure>
and observed wake length to predict the flexural rigidity values for each patch (Figure 13c) produces a set of values which have a similar range to those observed in the field. We note that the relationship in Figure 13a can produce multiple possible flexural rigidities for a given wake length. In this analysis, consistent with the majority of the field data (Figure 11c), the lower rigidity solution was used. The simplified model predicts a greater relative occurrence of higher rigidity patches compared to the field measurements (Figure 13c). This may reflect the dominance of higher rigidity shoots in determining patch wake length, above the influence of lower rigidity shoots also present within the patch. Given that the observed values in flexural rigidity vary over an order of magnitude, the simplified model's reproduction of this range does not provide direct evidence of the model's predictive capability but suggests the model may be capturing the key dynamics determining wake length.

# 4. Discussion

## 4.1. A Simplified Model for Vertical Shear Layers

In contrast to the CFD model, which provides detailed parameterization of shoot, leaf, and patch properties, as well as solution of the three-dimensional flow field, the simplified analytical model for vertical shear layers, based on equation (6), relies upon a single shear layer term, $\bar{U}/S_\delta\Delta U$ to encapsulate all of the flow and patch complexity. In this study, this parameterization was informed by the CFD simulations, but given approximate flow (velocity) and shoot (rigidity, diameter) data, site-specific models like that shown in Figure 13 could easily be produced following the same methodology and, except in the limiting case of very low velocity, would be expected to show the same overall trend with flexural rigidity.

While the model of Zong and Nepf (2012) produces a linear trend with increasing patch width, when considering the vertical shear layer, the relationship between wake length and shoot height is not linear, as the values of $x_d$ and $h_d$ are a function of flexural rigidity and determined by solving a beam equation. The model results (Figure 13a) can be understood by examining the theoretical model (Figure 1b) and the relative impacts of the terms $x_d$ and $h_d$ on wake length. It is clear from equation (6) that wake length increases with both $x_d$ and $h_d$; however, these two values are inversely related through the beam deflection model. The model results show a sharp increase in wake length at low rigidities, where the plant is in a more prone position, which is due to a large increase in $h_d$ as rigidity increases, compared to a small decrease in $x_d$ (Figure 13b). As flexural rigidity increases, the relative change in $h_d$ decreases and $x_d$ starts to decrease more significantly, reflecting the fact that in more upright positions, shoot reconfiguration leads to greater changes in $x_d$ than $h_d$ (Figure 13b). This causes a balancing out between the impacts of $h_d$ and $x_d$ and therefore the wake length reaches its maximum and then decreases slightly as the decrease in $x_d$ begins to dominate the increase in $h_d$ for rigidities above $EI= 3 \times 10^{-4} \text{ Nm}^2$. Due to the shear layer growth term, $S_\delta$, in equation (6), the influence of $h_d$ on $L_V$ is approximately 7 times that of $x_d$; hence, the reduction in wake length at higher rigidities is relatively small. In contrast, the CFD simulation data (black line in Figure 13a) demonstrates a larger decrease in wake length at higher rigidities, suggesting the influence of other factors which are not accounted for in the simplified analytical model. Such factors may include the additional impacts of streamlining on flow structure or changes to bleed flow and porosity within the patch as the plant reconfigures, but these factors seem to have most significant impact at higher rigidities. While the simplified model cannot capture such complexities, it is able to capture the impact of shoot reconfiguration and requires less than 0.0002% of the computational time compared to the CFD model enabling rapid investigation over a wide parameter space.

## 4.2. Dependence of Wake Length on the Vertical Shear Layer and Flexural Rigidity

The field data and CFD results on shear layer growth demonstrated that for pioneer saltmarsh vegetation (*Suaeda maritima*) patches, considering the ensemble-averaged patch geometry, patch height controls wake length rather than patch width. This geometric control is similar to what Liu et al. (2018) discovered for submerged rigid vegetation patches. There are, however, key differences with their study. First, their experiments were conducted with a high patch flow blockage, decreasing bleed flow and wake length (Chen et al., 2012; Liu et al., 2018). Second, here we consider flexible vegetation, and therefore as the plant reconfigures, $h_d$ and $x_d$ become significant. Finally, the variability in flexural rigidity across our data set confounds any direct scaling between patch height, $h$, and wake deposition length, $L_{TKE}$.
This is evident by looking at the individual field data for the observed vegetation patches. Assuming similar-
ity of shear layer growth as observed in the CFD simulations (Figure 10), it is possible to identify the theo-
retical controlling variable for each patch individually under rigid conditions. From the 45 patches, only
five exhibit theoretical dependence on patch width (w < 2h). We note that this classification does not account
for reconfiguration as this is not possible to quantify for the individual field patches. For the wakes deter-
mined by patch width (Figure 9b) there is a visual positive trend; however, there are too few points to pro-
duce a robust correlation. For the wakes determined by patch height (Figure 9c), there is less of an obvious
trend with wide scatter. This could be explained by variability in shoot height or density within the patch or
local velocity conditions. However, the simplified model demonstrates that it can also be explained by the
dependence of deflected patch height on flexural rigidity as shown in Figure 13c.

The field data for flexural rigidity show an order of magnitude of variation in values. Previous studies of both
freshwater and saltwater vegetation have also demonstrated large intraspecific variation in values for
flexural rigidity between samples on different stems from the same site (Feagin et al., 2011; Miler et al.,
2012; Paul et al., 2014; Rupprecht et al., 2015), and it is also likely that flexural rigidity varies depending
on patch exposure to the flow and the local age of vegetation patches (Anderson & Smith, 2014; Silinski
et al., 2018). The simplified model results (Figure 13c) demonstrate that the range of predicted flexural rigid-
ities agrees well with the observed range. However, they also demonstrate the large variation in wake length
caused by flexural rigidity which, over the range of flexural rigidity values observed, represents an uncer-
tainty in wake length of 1.55w. Considering the range of velocities that exhibited preferential deposition
(0.4–0.7 m/s), the simulation results suggest uncertainty in the wake length, due to local variations in velo-
city, of 1.1w (Figure 9a), which is also due in part to the impact of velocity on reconfiguration. By contrast,
the impact of shoot distribution on wake length was shown to produce a smaller uncertainty of 0.92w
(Figure 8). This highlights the relatively large impact of shoot reconfiguration on wake length, compared
to the within-patch structure. Furthermore, it demonstrates the importance of intraspecific variability in
flexural rigidity in determining wake length.

4.3. CFD-Biomechanical Model: Strengths, Limitations, and Future Directions

The coupled CFD-biomechanical model provides a new methodology for investigating patch reconfigura-
tion and flow-vegetation-sediment interaction, using high-resolution TLS and select field measurements
to parameterize a site-specific model. The model is not currently designed to represent the following:
(i) the complexity of stem and leaf scale flow-vegetation interactions, (ii) interpatch variation in geometric
and biomechanical properties, or (iii) specific patches within a site. Instead, it represents a site-averaged
vegetation patch model that enables consideration of complex shoot distribution and flexural rigidity
impacts which are either lacking in current models or restricted to single experimental cases (Tempest
et al., 2015).

Due to the number of simulations and associated computational time constraints, the simulations were run
as steady state RANS solutions. As such, while the model was able to capture static reconfiguration of the
vegetation to the time-averaged flow conditions, it could not capture the impacts of time-varying plant
motion. In reality, shoots undergo dynamic reconfiguration (Siniscalchi & Nikora, 2013) in response to
the turbulent flow which may act to either further dampen turbulence within the wake (Ghisalberti &
Nepf, 2006) or introduce higher levels of turbulent kinetic energy and turbulent events which may decrease
deposition (Ghisalberti & Nepf, 2009; Marjoribanks, Hardy, Lane, & Parsons et al., 2017).

Furthermore, Mont Saint-Michel is a dynamic environment which undergoes periodic tidal-driven sedimen-
tation events. Therefore, the observed sediment deposition is the product of multiple events and wake
deposition will in turn affect the flow field and therefore there may exist positive or negative feedbacks on
sediment deposition. Similarly, shoot-scale turbulence-induced resuspension and scour around the edges
of the patches may influence the wake flow and sediment processes (Tinoco & Coco, 2018; Yang & Nepf,
2018). Flow-sediment interactions will also drive, and be driven by, longer-term plant dynamics.
Individual plant life cycles, changes to patch structure, and colonization processes mean that patch density,
porosity, and flexural rigidity, three important variables within the model, will all vary over time. Moreover,
patches do not exist in isolation and therefore impacts of wave sheltering and larger-scale configuration will
affect flow velocities and lead to trait response.
The current CFD framework provides a mechanism for extending the capability of the model to capture these complexities through time-dependent solution using large eddy simulation, investigation of cumulative spatially distributed erosion and deposition over multiple tide cycles, and through multipatch simulation. Such developments would allow characterization of three-dimensional wake structure, permit investigation of the role of large-scale coherent turbulent structures on plant reconfiguration and sediment deposition, and enable investigation of wave damping and dissipation processes, which have been shown to be strongly dependent on flexural rigidity (Möller et al., 2014; Rupprecht et al., 2015). Nevertheless, the present model, with its representation of shoot reconfiguration and shoot-scale complexity, provides new insight into the impact of patch geometry, shoot distribution, and shoot biomechanics on wake length.

**5. Conclusions**

This study introduces a new field-parameterized vegetation patch model which permits the investigation of shoot-scale flow-vegetation interactions, bridging the gap between existing single-stem models and larger-scale rigid-cylinder or roughness approaches. The model was applied to investigate flow and potential sediment dynamics around pioneer saltmarsh vegetation patches. The main conclusions of this study are the following:

1. For the *Suaeda maritima* saltmarsh vegetation patches considered here, the vertical shear layer is dominant in controlling wake length, in contrast to much of the previous research on vegetation patches where width dominance is assumed.

2. Accordingly, flexural rigidity and patch reconfiguration become key controls on wake deposition length, in addition to patch height and width. Therefore, quantifying the intraspecific variability of biomechanical characteristics is critical to understanding the impact of vegetation on saltmarsh accretion.

3. Adapting the model of Zong and Nepf (2012) for height-controlled, flexible vegetation patches produces a simplified model for wake length which matches well with simulation data.

**Acknowledgments**

T. M. and R. J. H. acknowledge funding under NERC Grant NE/K003194/1, and D. L. and J. L. acknowledge project funding from CNRS/EC2CO/PNEC and PEPS terre-mer. Alix Mas and Lise Provost helped in field data collection. The authors would like to thank the Associate Editor and four anonymous reviewers whose comments have improved the paper. Data used to create the figures are available on the Loughborough University Data Repository (DOI:10.17028/rd.lboro.7588934).

**References**

Anderson, M. E., & Smith, J. M. (2014). Wave attenuation by flexible, idealized salt marsh vegetation. *Coastal Engineering, 83*, 82–92.
https://doi.org/10.1016/j.coastaleng.2013.10.004

Baptist, M. J., Babovic, V., Uthurburu, J. R., Keijzer, M., Uittenbogaard, R. E., Mynett, A., & Verwey, A. (2007). On inducing equations for vegetation resistance. *Journal of Hydraulic Research, 45*(4), 435–450. https://doi.org/10.1080/00221686.2007.9521778

Baranya, S., Olsen, N. R. B., Stoesser, T., & Sturm, T. (2012). Three-dimensional rans modeling of flow around circular piers using nested grids. *Engineering Applications of Computational Fluid Mechanics, 6*(4), 648–662. https://doi.org/10.1080/19942060.2012.11015449

Boothroyd, R. J., Hardy, R. J., Warburton, J., & Marjoribanks, T. I. (2017). Modeling complex flow structures and drag around a submerged plant of varied posture. *Water Resources Research, 53*, 2877–2901. https://doi.org/10.1002/2016WR020186

Brodu, N., & Lague, D. (2012). 3D terrestrial lidar data classification of complex natural scenes using a multi-scale dimensionality criterion: Applications in geomorphology. *ISPRS Journal of Photogrammetry and Remote Sensing, 68*(1), 121–134. https://doi.org/10.1016/j.isprsjprs.2012.01.006

Castellanos, E. M., Figueroa, M. E., & Davy, A. J. (1994). Nucleation and facilitation in saltmarsh succession: Interactions between *Spartina maritima* and *Arthrocnemum perenne*. *Journal of Ecology, 82*(2), 239–248. https://doi.org/10.2307/2261292

Chen, Z., Ortiz, A., Zong, L., & Nepf, H. (2012). The wake structure behind a porous obstruction and its implications for deposition near a finite patch of emergent vegetation. *Water Resources Research, 48*, W09517. https://doi.org/10.1029/2012WR012224

de Lima, P. H. S., Janzen, J. G., & Nepf, H. M. (2015). Flow patterns around two neighboring patches of emergent vegetation and possible implications for deposition and vegetation growth. *Environmental Fluid Mechanics, 15*(4), 881–898. https://doi.org/10.1007/s10652-015-9395-2

Dijkstra, J. T., & Uittenbogaard, R. E. (2010). Modeling the interaction between flow and highly flexible aquatic vegetation. *Water Resources Research, 46*, W12547. https://doi.org/10.1029/2010WR009246

Dimotakis, P. E. (1991). Turbulent free Shear Layer mixing and combustion. *High speed flight propulsion systems, 137*, 265–340. https://doi.org/10.2514/4.866104

Feagin, R. A., Irish, J. L., Möller, I., Williams, A. M., Colón-Rivera, R. J., & Mousavi, M. E. (2011). Short communication: Engineering properties of wetland plants with application to wave attenuation. *Coastal Engineering, 58*(3), 251–255. https://doi.org/10.1016/j.coastaleng.2010.10.003

General Public Licence (2014). Cloudcompare v2.5. General Public Licence. Retrieved from http://www.cloudcompare.org/

Ghisalberti, M., & Nepf, H. (2006). The structure of the shear layer in flows over rigid and flexible canopies. *Environmental Fluid Mechanics, 6*(3), 277–301. https://doi.org/10.1007/s10652-006-0002-4

Ghisalberti, M., & Nepf, H. (2009). Shallow flows over a permeable medium: The hydrodynamics of submerged aquatic canopies. *Transport in Porous Media, 78*(3 SPEC. ISS), 385–402. https://doi.org/10.1007/s11242-009-9434-x

Hamann, E., & Puijalon, S. (2013). Biomechanical responses of aquatic plants to aerial conditions. *Annals of Botany, 112*(9), 1869–1878.
https://doi.org/10.1093/aob/mct221

Hamed, A. M., Sadowski, M. J., Nepf, H. M., & Chamorro, L. P. (2017). Impact of height heterogeneity on canopy turbulence. *Journal of Fluid Mechanics, 813*, 1176–1196. https://doi.org/10.1017/jfm.2017.22
Hardy, R. J., Lane, S. N., Lawless, M. R., Best, J. L., Elliott, L., & Ingham, D. B. (2005). Development and testing of a numerical code for treatment of complex river channel topography in three-dimensional CFD models with structured grids. *Journal of Hydraulic Research*, 43(5), 468–480. https://doi.org/10.1080/00221680509500145

Hu, Z., Lei, J., Liu, C., & Nepf, H. (2018). Wake structure and sediment deposition behind models of submerged vegetation with and without flexible leaves. *Advances in Water Resources*, 118, 28–38. https://doi.org/10.1016/j.adwatres.2018.06.001

Lague, D., Brodu, N., & Leroux, J. (2013). Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (N-Z). *ISPRS Journal of Photogrammetry and Remote Sensing*, 82, 10–26. https://doi.org/10.1016/j.isprsjprs.2013.04.009

Langlois, E., Bonis, A., & Bouzillé, J. B. (2003). Sediment and plant dynamics in saltmarshes pioneer zone: *Puccinellia maritima* as a key species? *Estuarine, Coastal and Shelf Science*, 56(2), 239–249. https://doi.org/10.1016/S0272-7714(02)00185-3

Leroux, J. (2013). Chenaux tidaux et dynamique des prés-salés en régime méga-tidal: Approche multi-temporelle du siècle à l'événement de marée. Université Rennes 1; Université Européenne de Bretagne.

Liu, C., Hu, Z., Lei, J., & Nepf, H. (2018). Vortex structure and sediment deposition in the wake behind a finite patch of model submerged vegetation. *Journal of Hydraulic Engineering*, 144(2), 04017065. https://doi.org/10.1061/(ASCE)HY.1943-7900.0001408

Luhar, M., & Nepf, H. M. (2011). Flow-induced reconfiguration of buoyant and flexible aquatic vegetation. *Limnology and Oceanography*, 56(6), 2003–2017. https://doi.org/10.4319/lo.1011.56.6.2003

Marjoribanks, T. I., Hardy, R. J., Lane, S. N., & Parsons, D. R. (2014). High-resolution numerical modelling of flow-vegetation interactions. *Journal of Hydraulic Research*, 52(6), 775–793. https://doi.org/10.1080/00221686.2014.948502

Marjoribanks, T. I., Hardy, R. J., Lane, S. N., & Parsons, D. R. (2017). Does the canopy mixing layer model apply to highly flexible aquatic vegetation? *Insights from numerical modelling*. *Environmental Fluid Mechanics*, 17(2), 277–301. https://doi.org/10.1007/s10652-016-9482-z

Marjoribanks, T. I., Hardy, R. J., Lane, S. N., & Tancock, M. J. (2017). Patch-scale representation of vegetation within hydraulic models. *Earth Surface Processes and Landforms*, 42(5), 699–710. https://doi.org/10.1002/esp.4015

Miler, O., Albayrak, I., Nikora, V., & O'Hare, M. (2012). Biomechanical properties of aquatic plants and their effects on plant-flow interactions in streams and rivers. *Aquatic Sciences*, 74(1), 31–44. https://doi.org/10.1007/s00027-011-0188-5

Möller, I. (2006). Quantifying saltmarsh vegetation and its effect on wave height dissipation: Results from a UK East coast saltmarsh. *Estuarine, Coastal and Shelf Science*, 69(3–4), 337–351. https://doi.org/10.1016/j.ecss.2006.05.003

Möller, I., Kudella, M., Rupprecht, F., Spencer, T., Paul, M., van Wesenbeeck, B. K., et al. (2014). Wave attenuation over coastal salt marshes under storm surge conditions. *Nature Geoscience*, 7(10), 727–731. https://doi.org/10.1038/NGEO2251

Nepf, H. M. (1999). Drag, turbulence, and diffusion in flow through emergent vegetation. *Water Resources Research*, 35(2), 479–489. https://doi.org/10.1029/1998WR900069

Nepf, H. M., & Vivoni, E. R. (2000). Flow structure in depth-limited, vegetated flow. *Journal of Geophysical Research*, 105(C12), 28,547–28,557. https://doi.org/10.1029/2000JC900145

Neumeier, U., & Amos, C. L. (2006). The influence of vegetation on turbulence and flow velocities in European salt-marshes. *Sedimentology*, 53(2), 259–277. https://doi.org/10.1111/j.1365-3091.2006.00772.x

Nicolle, A., & Eames, I. (2011). Numerical study of flow through and around a circular array of cylinders. *Journal of Fluid Mechanics*, 679, 1–31. https://doi.org/10.1017/jfm.2011.77

Niklas, K. J. (1992). *Plant biomechanics: An engineering approach to plant form and function*. Chicago: University of Chicago Press.

Panton, R. L. (1996). *Incompressible flow*, (2nd ed.). New York: John Wiley and Sons.

Paul, M., Henry, P. Y. T., & Thomas, R. E. (2014). Geometrical and mechanical properties of four species of northern European brown macroalgae. *Coastal Engineering*, 84, 73–80. https://doi.org/10.1016/j.coastaleng.2013.11.007

Pope, N. D., Widdows, J., & Brinsley, M. D. (2006). Estimation of bed shear stress using the turbulent kinetic energy approach—A comparison of annular flume and field data. *Continental Shelf Research*, 26(8), 959–970. https://doi.org/10.1016/j.csr.2006.02.010

Pope, S. B. (2000). Turbulent flows. *Book*, 1. https://doi.org/10.1088/1468-5248/1/1/702

Rodi, W. (1980). Turbulence models and their application in hydraulics—A state of the art review. *IAHR Publication*. Retrieved from http://adsabs.harvard.edu/abs/1980STIA...8121395R

Rupprecht, F., Möller, I., Evans, B., Spencer, T., & Jensen, K. (2015). Biophysical properties of salt marsh canopies—Quantifying plant stem flexibility and above ground biomass. *Coastal Engineering*, 100, 48–57. https://doi.org/10.1016/j.coastaleng.2015.03.009

Salehi, M., & Strom, K. (2011). Using velocimeter signal to noise ratio as a surrogate measure of suspended mud concentration. *Continental Shelf Research*, 31(9), 1020–1032. https://doi.org/10.1016/j.csr.2011.03.008

Shepard, C. C., Crain, C. M., & Beck, M. W. (2011). The protective role of coastal marshes: A systematic review and meta-analysis. *PLoS ONE*, 6(11). https://doi.org/10.1371/journal.pone.0027374

Shi, Y., Jiang, B., & Nepf, H. M. (2016). Influence of particle size and density, and channel velocity on the deposition patterns around a circular patch of model emergent vegetation. *Water Resources Research*, 52, 1044–1055. https://doi.org/10.1002/2015WR018278

Shields, A. (1936). Anwendung der Aehnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung. Mitt. Preuss. Versuchsanst. Wasserbau Schiffbau. https://doi.org/10.1061/(ASCE)0733-9429(1995)121:11/766

Silinski, A., Schoutens, K., Puijalon, S., Schoelynck, J., Luyckx, D., Troch, P., et al. (2018). Coping with waves: Plasticity in tidal marsh plants as self-adapting coastal ecosystem engineers. *Limnology and Oceanography*, 63(2), 799–815. https://doi.org/10.1002/lno.10671

Siniscalchi, F., & Nikora, V. (2013). Dynamic reconfiguration of aquatic plants and its interrelations with upstream turbulence and drag forces. *Journal of Hydraulic Research*, 51(1), 46–55. https://doi.org/10.1080/00221686.2012.743486

Soulsby, R. (1997). *Dynamics of marine sands: A manual for practical applications*. London: Thomas Telford. https://doi.org/10.1680/ doms.25844

Soulsby, R. L. (1981). Measurements of the reynolds stress components close to a marine sand bank. *Marine Geology*, 42(1-4), 35–47. https://doi.org/10.1016/S0707-4571(08)70293-0

Spalding, D. B. (1980). Mathematical modelling of fluid mechanics, heat transfer and mass transfer processes. London.

Sukhodolova, T. A., & Sukhodolov, A. N. (2012). Vegetated mixing layer around a finite-size patch of submerged plants: 1. Theory and field experiments. *Water Resources Research*, 48, W12506. https://doi.org/10.1029/2011WR011804

Temmerman, S., Bouma, T. J., Govers, G., Wang, Z. B., De Vries, M. B., & Herman, P. M. J. (2005). Impact of vegetation on flow routing and sedimentation patterns: Three-dimensional modeling for a tidal marsh. *Journal of Geophysical Research*, 110, F04019. https://doi.org/10.1029/2005JF000301

Temmerman, S., Meire, P., Bouma, T. J., Herman, P. M. J., Ysebaert, T., & De Vriend, H. J. (2013). Ecosystem-based coastal defence in the face of global change. *Nature*, 504(7478), 79–83. https://doi.org/10.1038/nature12859
Tempest, J. A., Möller, I., & Spencer, T. (2015). A review of plant-flow interactions on salt marshes: The importance of vegetation structure and plant mechanical characteristics. *Wiley Interdisciplinary Reviews: Water, 2*(6), 669-681. https://doi.org/10.1002/wat2.1103

Tessier, M., Gloaguen, J. C., & Lefeuvre, J. C. (2000). Factors affecting the population dynamics of *Suaeda maritima* at initial stages of development. *Plant Ecology*, **147**(2), 193-203. https://doi.org/10.1023/A:1009854204841

Tinoco, R. O., & Coco, G. (2018). Turbulence as the main driver of resuspension in oscillatory flow through vegetation. *Journal of Geophysical Research: Earth Surface*, **123**, 891–904. https://doi.org/10.1002/2017JF004504

Tison, J. M., & de Foucault, B. (2014). Flora Gallica: Flore de France. Biotope.

Wu, G., Li, H., Liang, B., Shi, F., Kirby, J. T., & Mieras, R. (2017). Subgrid modeling of salt marsh hydrodynamics with effects of vegetation and vegetation zonation. *Earth Surface Processes and Landforms*, **42**(12), 1755–1768. https://doi.org/10.1002/esp.4121

Yang, J. Q., Chung, H., & Nepf, H. M. (2016). The onset of sediment transport in vegetated channels predicted by turbulent kinetic energy. *Geophysical Research Letters*, **43**, 11,261–11,268. https://doi.org/10.1002/2016GL071092

Yang, J. Q., & Nepf, H. M. (2018). A turbulence-based bed-load transport model for bare and vegetated channels. *Geophysical Research Letters*, **45**, 10,428–10,436. https://doi.org/10.1029/2018GL079319

Zong, L., & Nepf, H. (2012). Vortex development behind a finite porous obstruction in a channel. *Journal of Fluid Mechanics*, **691**, 368–391. https://doi.org/10.1017/jfm.2011.479