1. Introduction
Wall-mounted flexible structures, such as seagrasses, mangroves, kelps, forests and coral reefs, are ubiquitous in natural aquatic environments, and play crucial roles in ecosystem functioning (Carpenter & Lodge Reference Carpenter and Lodge1986; Green & Short Reference Green and Short2003; Unsworth et al. Reference Unsworth, Cullen-Unsworth, Jones and Lilley2022). Aquatic vegetation provides habitats and regulates oxygen, carbon and nutrient cycling (Wilcock et al. Reference Wilcock, Champion, Nagels and Croker1999; Schulz et al. Reference Schulz, Kozerski, Pluntke and Rinke2003). Within vegetated meadows, nutrients and sediments tend to be retained and deposited, potentially improving water quality, supplying food sources for aquatic organisms, and contributing to bed stabilisation (Carpenter & Lodge Reference Carpenter and Lodge1986; Sand-Jensen Reference Sand-Jensen1998). From an engineering and environmental perspective, submerged canopies are also effective in attenuating wave energy and modulating sediment transport, contributing to coastal stability (Day et al. Reference Day2007; Temmerman et al. Reference Temmerman, Meire, Bouma, Herman, Ysebaert and De Vriend2013). However, despite their ecological and hydrodynamic importance, the coupled dynamics between vegetation motion and surrounding flow structures remains insufficiently understood. The underlying physical mechanisms depend strongly on vegetation properties and spatial arrangements (Jackson & Winant Reference Jackson and Winant1983; Finnigan Reference Finnigan2000; Rominger & Nepf Reference Rominger and Nepf2011). Accurate characterisation of flow–structure interactions (FSIs) within submerged canopies requires spatially and temporally resolved measurements of the in-canopy flow, vegetation motion and associated hydrodynamic forces, which are inherently difficult to obtain experimentally. High-fidelity numerical simulation therefore serves as a powerful tool to resolve the coupled fluid and structural dynamics.
The flow characteristics within and above rigid vegetation canopies are governed by their geometric configuration and structural properties (Nepf Reference Nepf, Ghisalberti, White and Murphy2012). The canopy geometry is characterised by the solid volume fraction (
$\phi$
) occupied by canopy elements. When the element is a circular cylinder,
$\phi\approx0.25\pi d^2/D^2$
, where
$d$
and
$D$
denote the vegetation diameter and array width, respectively (Nepf Reference Nepf, Ghisalberti, White and Murphy2012; Nicholas et al. Reference Nicholas, Omidyeganeh and Pinelli2022). Based on
$\phi$
, the value of which varies with vegetation morphology, canopies are categorised as sparse, dense and transitional (Belcher et al. Reference Belcher, Jerram and Hunt2003; Coceal & Belcher Reference Coceal and Belcher2004; Nepf Reference Nepf, Ghisalberti, White and Murphy2012). Sparse canopies resemble rough-wall boundary layers without velocity inflection (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004; Nepf et al. Reference Nepf2007), whereas dense canopies develop two inflection points associated with enhanced shear layers (Nepf Reference Nepf, Ghisalberti, White and Murphy2012). Transitional canopies exhibit a single inflection point that separates the flow into inner- and outer-layer regions. The inner region (
$z/h\ll1$
) is dominated by a von Kármán vortex street, while the outer layer (
$z/h\gt2$
) behaves like a rough-wall boundary layer. Between them, a turbulent mixing layer develops at the inflection point, where Kelvin–Helmholtz (K–H) instabilities generate large-scale spanwise vortices (Nicholas et al. Reference Nicholas, Omidyeganeh and Pinelli2022).
Flexible vegetation exhibits additional complexity due to FSI (Sriram et al. Reference Sriram, Stoesser, Yan and Murali2023). The dynamic responses depend strongly on flexibility, density and morphology (Kim et al. Reference Kim, Yu and Kaplan2024). Flexibility is commonly expressed as
$K = \textit{EI}/ 0.5 \rho d h^3 U_{\textit{ref}}^2$
, where
$E$
is Young’s modulus,
$I$
is the cross-sectional momentum of inertia,
$\rho$
is the fluid density,
$d$
is the side length or diameter, and
$U_{\textit{ref}}$
is the reference velocity. Flexibility affects flow instability in two competing ways (Wong et al. Reference Wong, Trinh and Chapman2020): (i) reconfiguration of vegetation reduces drag and stabilises the flow, whereas (ii) the decoupling of fluid and structural perturbations may enhance instability. Depending on flexibility, canopy motion ranges from rigid to strongly prone state, with ‘monami’ representing coherent waving induced by K–H vortices (Ackerman & Okubo Reference Ackerman and Okubo1993; Carollo et al. Reference Carollo, Ferro and Termini2005; Okamoto et al. Reference Okamoto, Nezu and Sanjou2016). Subsequent studies (Tschisgale et al. Reference Tschisgale, Löhrer, Meller and Fröhlich2021; Wang et al. Reference Wang, Cheng and Ting2022a
; Monti et al. Reference Monti, Olivieri and Rosti2023; Ni et al. Reference Ni, Ji, Xu, Zhang and Liang2023; Ni et al. Reference Ni, Zhang, Ji, Xu and Zhang2025) demonstrated that monami persists across a range of flexibilities and primarily arises from large-scale turbulent fluctuations near the canopy top.
Flow within vegetated canopies contains both canopy-scale and stem-scale vortices. The shear and velocity inflection at the canopy top trigger K–H instabilities that generate large coherent vortices (Raupach et al. Reference Raupach, Finnigan and Brunet1996; Ghisalberti & Nepf Reference Ghisalberti and Nepf2002; Singh et al. Reference Singh, Bandi, Mahadevan and Mandre2016; Wong et al. Reference Wong, Trinh and Chapman2020; Vieira et al. Reference Vieira, Allshouse and Mahadevan2023). These vortices enhance turbulent mixing and vertical momentum transport, promoting the intrusion of Reynolds shear stress into the canopy interior, and modifying the in-canopy flow dynamics (Ghisalberti & Nepf Reference Ghisalberti and Nepf2006; Wu et al. Reference Wu, Wen, Yen, Weng and Wang2023). Stem-scale vortices, by contrast, form in the wake of individual elements (Nepf & Vivoni Reference Nepf and Vivoni2000; Ghisalberti & Nepf Reference Ghisalberti and Nepf2009; Nepf Reference Nepf, Ghisalberti, White and Murphy2012; He et al. Reference He, Liu and Shen2022).
The velocity distribution within canopies is also influenced by vegetation arrangement, height and density. Flow development proceeds from an adjustment region, where the velocity decays exponentially from the canopy edge (Chen et al. Reference Chen and Li2013), to a fully developed regime characterised by equilibrium between drag and turbulence production (Belcher et al. Reference Belcher, Jerram and Hunt2003; Rominger & Nepf Reference Rominger and Nepf2011; Lei & Nepf Reference Lei and Nepf2021; Zhang et al. Reference Zhang, Wang, Li, Guo, Li, Xu, Yu and Wang2024). While models for rigid canopies are well validated, the modulation of velocity and turbulence by stem flexibility, particularly for highly deformable vegetation, remains insufficiently understood. This gap motivates the present study, which aims to elucidate how vegetation flexibility reshapes flow structures and hydrodynamic forces within and above canopy arrays.
Hydrodynamic forces govern flow resistance, energy dissipation and sediment transport in vegetated flows. The vegetation-induced drag coefficient is commonly defined as
$C_d = F_d / 0.5 \rho U^2 {d} h$
, where
$F_d$
is the drag force acting on the vegetation, and
$U$
is the flow velocity. Early resistance models rely on empirical bed-stress or Manning-type formulations (Thompson & Roberson Reference Thompson and Roberson1976; Arcement et al. Reference Arcement, Schneider and Verne1988; Kadlec Reference Kadlec1990). These models assume energy loss solely through bed stress, and require empirical calibration, limiting their accuracy for vegetated flows. More recent approaches employ a bulk drag coefficient based on bulk or pore velocity (Wu et al. Reference Wu, Wu, Wu, Zhang, Feng, Zhang and Zhao1999; Zhao et al. Reference Zhao and Nepf2013). However, these reference velocities do not fully represent the effective flow acting on the canopy elements because canopy drag depends on resolved FSI; purely experimental estimates remain indirect (Luhar et al. Reference Luhar, Coutu, Infantes, Fox and Nepf2010; Pujol et al. Reference Pujol, Serra, Colomer and Casamitjana2013).
The prediction of sediment transport in vegetated flows is equally critical. Traditional formulations (Shields Reference Shields1936; Einstein Reference Einstein1942; Meyer-Peter & Müller Reference Meyer-Peter and Müller1948) relate bedload transport to the mean bed shear stress, whereas recent studies have shown that turbulence generated by vegetation substantially enhances sediment mobility (Heathershaw & Thorne Reference Heathershaw and Thorne1985; Sumer et al. Reference Sumer, Chua, Cheng and Fredsøe2003; Diplas et al. Reference Diplas, Dancey, Celik, Valyrakis, Greer and Akar2008; Shih & Diplas Reference Shih, Wang and Ting2018; Tseng & Tinoco Reference Tseng and Tinoco2021). In vegetated channels, bare-bed models for bed shear stress (
$\tau _b$
) often produce substantial errors, since bed stress alone does not fully account for vegetation generated turbulence (Zhao & Nepf Reference Zhao, Cheng, Wang and Tan2021; Yang et al. Reference Yang, Zhang, Qi and Li2025). Sediment transport can be underestimated by several orders of magnitude when vegetation-induced turbulence is neglected (Yang & Nepf Reference Yang and Nepf2018). Both the near-bed shear stress and the turbulent kinetic energy (TKE) should be considered in sediment transport models (Tinoco & Cowen Reference Tinoco and Cowen2013; Tinoco & Coco Reference Tinoco and Coco2014; Tinoco & Coco Reference Tinoco and Coco2016; Tinoco & Coco Reference Tinoco and Coco2018; Armanini & Cavedon Reference Armanini and Cavedon2019; Zhao & Nepf Reference Zhao, Cheng, Wang and Tan2021; Wang et al. Reference Wang, He, Dey and Fang2023). Numerical simulations additionally enable direct evaluation of both the mean and root mean square (r.m.s.) bed shear stress, providing a quantitative measure of sediment entrainment (Koken & Constantinescu Reference Koken and Constantinescu2021).
Within emergent canopies, turbulence is primarily generated by wakes of individual stems and branches (Nepf Reference Nepf, Ghisalberti, White and Murphy2012). The magnitude and distribution of this vegetation-generated turbulence depend on stem density, local Reynolds number and plant morphology (Yager & Schmeeckle Reference Yager and Schmeeckle2013; Shan et al. Reference Shan, Zhao, Liu and Nepf2020; Xu & Nepf Reference Xu and Zhou2020; Zhao & Nepf Reference Zhao, Cheng, Wang and Tan2021). In most cases, vegetation-generated turbulence dominates bed-generated shear stress, making it the primary driver of sediment motion (Yang et al. Reference Yang, Chung and Nepf2016; Yang & Nepf Reference Yang and Nepf2018; Yang & Nepf Reference Yang and Nepf2019). For submerged canopies, however, an additional source of turbulence arises from the drag discontinuity at the canopy top, which induces a stronger shear layer and promotes turbulence production (King et al. Reference King, Tinoco and Cowen2012; Zhang et al. Reference Zhang, Lei, Huai and Nepf2020). Lateral shear layers are also important because vortices generated near the canopy edges can penetrate towards the bed. As turbulence modulates sediment mobility and energy dissipation, quantifying the role of canopy-generated turbulence remains of both theoretical and practical significance.
Stem reconfiguration alters the effective frontal area, redistributes drag, and modifies turbulence and mixing. Understanding these coupled processes is therefore essential, and forms a central objective of the present study. Despite considerable previous research, simultaneously resolving all relevant quantities, including stem dynamics, wake evolution, flow distribution, sediment transport and energy dissipation, at high spatial and temporal resolution remains challenging. In this work, we conduct direct numerical simulations of flow through a square array of flexible aquatic stems, representing an idealised model of submerged vegetation with uniform diameter and flexibility along the vertical axis. A newly developed numerical framework, combining the immersed boundary method with the vector form intrinsic finite element (VFIFE) method (Chen & Li Reference Chen, Xu and Ji2025), is employed. We investigate how vegetation flexibility modulates FSI and hydrodynamics, with objectives to: (i) characterise flow structures above and within flexible arrays; (ii) quantify stem structural responses under hydrodynamic loading; (iii) assess the effects of flexibility on hydrodynamic drag and in-canopy flow redistribution, including bed shear stress and TKE; and (iv) elucidate the physical mechanisms governing structural dynamics and hydrodynamic loading.
The remainder of this paper is organised as follows. Section 2 introduces the numerical methodology and simulation set-up. Section 3 presents the stem dynamics, three-dimensional (3-D) vortex structures, velocity distribution, hydrodynamic forces, bed shear stress and TKE variations with stem flexibility. The underlying mechanisms governing the stem motion and hydrodynamic loading are also discussed in § 4. Section 5 summarises the main findings and implications for coastal protection and sediment transport.
2. Numerical methodology
2.1. Numerical method
The numerical simulations in this study were conducted using an in-house computational fluid dynamics code. The governing equations for the incompressible flow are the Navier–Stokes and continuity equations expressed as
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}&=-\frac {1}{\rho } \boldsymbol{\nabla }\!p+\nu\, {\nabla} ^2 \boldsymbol{u}+\boldsymbol{f},\\ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0,\nonumber \end{align}
where
$\boldsymbol{u}$
is the velocity of the fluid flow,
$t$
is time,
$p$
is pressure,
$\rho$
is the fluid density,
$\nu$
is the kinematic viscosity, and
$\boldsymbol{f}$
is the extra body force that represents the action of the immersed objects on the fluid.
To discretise (2.1), the second-order Adams–Bashforth time scheme and the second-order central difference scheme are applied, which lead to the conservation form
\begin{align} &\boldsymbol{u}^{n+1}=\boldsymbol{u}^n+\delta t\left (\frac {3}{2} \boldsymbol{h}^n-\frac {1}{2} \boldsymbol{h}^{n-1}-\frac {3}{2} \boldsymbol{\nabla }\!p^n+\frac {1}{2} \boldsymbol{\nabla }\!p^{n-1}\right )+\boldsymbol{f}^{n+{1}/{2}} \delta t, \\ &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{n+1}=0,\nonumber \end{align}
where
$\boldsymbol{h}=\boldsymbol{\nabla }\boldsymbol{\cdot }(-\boldsymbol{u} \boldsymbol{u}+\nu (\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^{\rm T} ) )$
is composed of the convective and diffusive terms, in which the superscript
${\rm T}$
indicates the transpose of a matrix. The superscript
$n$
indicates the time step. The extra force
$\boldsymbol{f}$
is determined by enforcing the velocity-consistency condition
$\boldsymbol{U}^{n+1} = \boldsymbol{V}^{n+1}$
on the fluid–structure interface, where
$\boldsymbol{U}^{n+1}$
is the interpolated fluid velocity, and
$\boldsymbol{V}^{n+1}$
is the update structure velocity. The extra force
$\boldsymbol{f}$
is formally evaluated at the time step
$n+1/2$
:
\begin{equation} \boldsymbol{f}^{n+{1}/{2}} \delta t=D (\boldsymbol{F}^{n+{1}/{2}} )=D\left (\boldsymbol{V}^{n+1}-I\left (u^n+\delta t\left (\frac {3}{2} \boldsymbol{h}^n-\frac {1}{2} \boldsymbol{h}^{n-1}-\frac {3}{2} \boldsymbol{\nabla }\!p^n+\frac {1}{2} \boldsymbol{\nabla }\!p^{n-1}\right )\right )\right )\!, \end{equation}
where
$\boldsymbol{F}$
is the extra body force on the structure. The interpolation and distribution functions
$I$
and
$D$
, suggested in Peskin (Reference Peskin1972), are applied, and are responsible for data transfer between the non-conforming Eulerian and Lagrangian grids.
The dynamics of a wall-mounted flexible stem is simulated using the VFIFE method (Shih et al. Reference Shih and Diplas2004; Ting et al. Reference Ting, Shih and Wang2004a ,Reference Ting, Shih and Wang b , Reference Ting, Duan and Wu2012; Ting & Wang Reference Ting and Wang2008). The VFIFE method is an analytical-numerical approach in structural mechanics based on point value representation and vector form theory. By employing a lumped mass matrix, the method ensures that external forces applied at one spatial point do not affect the motion of other points, thereby eliminating the need for constructing stiffness matrices or solving large nonlinear systems. Compared with conventional finite element methods, the VFIFE method effectively handles large deformations, collisions and elastoplastic behaviours, while avoiding common nonlinear issues such as non-convergence and matrix singularity. Based on the point value description, the stem is discretised into a series of spatial nodes connected by beam elements. The motion and deformation of the stem are represented by the displacement of these nodes, each of which follows Newton’s second law of motion:
\begin{align} \begin{aligned} &m \dot {\boldsymbol{V}}+\zeta m \boldsymbol{V}+\boldsymbol{F}^{\textit{int}}=\boldsymbol{F}^{\textit{ext}}, \\ &m \dot {\boldsymbol{\omega }}+\zeta m \boldsymbol{\omega }+\boldsymbol{M}^{\textit{int}}=\boldsymbol{M}^{\textit{ext}}, \end{aligned} \end{align}
where
$m$
and
$I$
are the mass and momentum of inertia matrix of the node, respectively. Here,
$\dot {\boldsymbol{V}}$
,
$\boldsymbol{V}$
,
$\dot {\boldsymbol{\omega }}$
and
$\boldsymbol{\omega }$
are the acceleration, velocity, angular acceleration and angular velocity of the node, respectively,
$\boldsymbol{F}^{\textit{int}}$
,
$\boldsymbol{F}^{\textit{ext}}$
,
$\boldsymbol{M}^{\textit{int}}$
and
$\boldsymbol{M}^{\textit{ext}}$
are the internal force, external force, internal moment and external moment, respectively, and
$\zeta$
is the structural damping coefficient, which is zero in the present study. An explicit iterative scheme based on the central difference is adopted to solve (2.4) as follows:
\begin{align} \begin{aligned} &\boldsymbol{x}^{n+1}=\Delta t_s^2 \frac {\boldsymbol{F}^{\textit{ext}}-\boldsymbol{F}^{\textit{int}}}{m}+2 \boldsymbol{x}^n-\boldsymbol{x}^{n-1},\\ &\boldsymbol{\alpha }^{n+1}=\Delta t_s^2 \frac {\boldsymbol{M}^{\textit{ext }}-\boldsymbol{M}^{\textit{int }}}{m}+2 \boldsymbol{\alpha }^n-\boldsymbol{\alpha }^{n-1},\\ &\boldsymbol{V}^{n+1}=\frac {\boldsymbol{x}^{n+1}-\boldsymbol{x}^n}{\Delta t_s}, \quad \omega ^{n+1}=\frac {\boldsymbol{\alpha }^{n+1}-\boldsymbol{\alpha }^n}{\Delta t_s}, \end{aligned} \end{align}
where
$\boldsymbol{x}$
,
$\boldsymbol{\alpha }$
and
$t_s$
represent the node position, rotation angle of the node, and the structural time step, respectively. The temporal trajectory of any node on the stem is described by a set of spatial points, and the process between the spatial points is defined as the path element. Proper division of the path element allows the properties of each node to vary between adjacent elements. Collisions between neighbouring vegetation or between the tip and the wall do not occur, therefore no contact model is employed in the simulations.
(
$a$
) Computational domain and (
$b$
) configuration of the flexible stem array. The array consists of 25 stems arranged in a square array with uniform spacing
$4d$
in both the streamwise and transverse directions. Each stem has height-to-diameter ratio
$h/d = 10$
.
2.2. Simulation set-up for flow through flexible stem arrays
In the present simulations, a square patch comprising 25 flexible stems is modelled, as shown in figure 1(
$a$
), representing a short canopy patch (Lei & Nepf Reference Lei and Nepf2021). Flow within such a patch remains in the adjustment region, where the streamwise velocity decays from the canopy edge and is generally higher than that in a fully developed region over longer canopies (Chen et al. Reference Chen and Li2013). Each stem has a height of
$10d$
, and the stems are uniformly spaced
$4d$
apart in both the streamwise and transverse directions, as shown in figure 1(
$b$
). The stem is discretised into 160 beam elements along its height, resulting in 161 circular slices. Each slice is represented by 256 uniformly distributed immersed boundary points (IBPs), ensuring at least one IBP per grid cell, and maintaining adequate spatial resolution. Accordingly, the vertical grid resolution is
$\Delta z / d = 1/16$
.
The computational domain extends
$40d$
,
$30d$
and
$16d$
in the streamwise (
$x$
), transverse (
$y$
) and vertical (
$z$
) directions, respectively. The first stem row is located
$8d$
downstream of the inlet, and the lateral clearance between the outer stems and the side boundaries is
$7d$
, as shown in figure 1(
$b$
). The upstream distance is selected to allow the imposed inlet condition to adjust before reaching the array. A shorter inlet length may lead to slight variations in the approaching flow. The downstream distance primarily affects the development of the wake and coherent vortical structures near the outlet, whereas its influence on the mean hydrodynamic forces acting on the stems is negligible. Previous domain-independence studies showed that extending the outlet further changes the mean forces by less than
$3\,\%$
(Chen et al. Reference Chen, Ji, Alam, Xu, An, Tong and Zhao2015; Chen et al. Reference Chen and Li2022; Tu et al. Reference Tu, Yeoh and Liu2018). The submergence ratio (water depth/stem height) is 1.6, corresponding to shallow submergence (Nepf Reference Nepf, Ghisalberti, White and Murphy2012). To cover the potential motion of the stems, a high-resolution uniform mesh covers
$32d\times26d$
in the central region (
$3d\leqslant x\leqslant35d$
,
$2d\leqslant y\leqslant28d$
), with
$\Delta x / d=\Delta y / d=1/32$
, as indicated by the blue rectangle in figure 1(
$b$
). Outside this zone, a stretched grid is applied to reduce computational costs. The mesh resolution follows the criteria validated in our previous work (Chen & Li Reference Chen, Xu and Ji2025). Boundary conditions consist of a Dirichlet inflow (
$u = U_{\infty }$
,
$v = 0$
,
$w= 0$
), a Neumann outflow (
$\partial u / \partial x = 0$
,
$\partial v / \partial x = 0$
,
$\partial w / \partial x = 0$
), free-slip lateral boundaries (
$\partial u / \partial y = 0$
,
$ v = 0$
,
$\partial w / \partial y = 0$
), and a no-slip bottom wall (
$u = 0$
,
$v = 0$
,
$w= 0$
). The free-slip treatment of the upper boundary (
$\partial u / \partial z = 0$
,
$\partial v / \partial z = 0$
,
$ w = 0$
) neglects free-surface effects. This approximation is justified because the vortical structures above the canopy are confined to a limited vertical extent (typically within
$3d$
for cases without K–H vortices), which is substantially smaller than the
$6d$
clearance between the canopy top and the upper boundary adopted in the present simulations. Saha (Reference Saha2013) demonstrated that for flow over finite wall-mounted cylinders, a clearance exceeding
$5d$
is sufficient to eliminate noticeable free-surface effects on the flow field. In addition, stem reconfiguration further reduces the effective canopy height, increasing the distance from the free surface. Therefore, free-surface effects are expected to be negligible and are not explicitly resolved in the present simulations.
The flexibility of aquatic vegetation in nature typically ranges from
$K=2.0 \times 10^{-4}$
to
$K=3.0 \times 10^{2}$
(Luhar & Nepf Reference Luhar and Nepf2011). Our previous study (Chen & Li Reference Chen, Xu and Ji2025) showed that stems with
$K \gt 1.0$
behave as rigid bodies; hence only cases with
$K \leqslant1.0$
are considered here. To systematically investigate the effect of flexibility on the flow and stem dynamics, a broad range of
$K$
is considered, with cases corresponding to
$K = 1$
, 0.32, 0.1, 0.032, 0.01 and 0.0032. The stem density is set to
$\rho _{s} / \rho = 2$
. The simulations were performed on the computational facilities of the National Supercomputing Centre, Singapore. The non-dimensional time step is set to
$2 \times 10^{-5}$
, yielding a CFL number of approximately
$10^{-3}$
. Each case was computed in parallel using 256 CPU cores, and required approximately 3–4 months of wall-clock time per case.
3. Results
This section presents the stem-array dynamics, stem deformation, vortex dynamics, inter-stem interactions, flow distributions, hydrodynamic loading and bed shear stress as functions of stem flexibility. Table 1 summarises the key parameters of the investigated cases. Three distinct response regimes are identified within the simulated range of
$K$
. For brevity, only the representative cases
$K = 1$
, 0.32, 0.1 and 0.032 are presented, as they capture all three response regimes.
Main geometric and flow parameters of the investigated cases. The array width
$D$
is defined as the distance between the outer edges of the outermost stems, i.e.
$D = (N-1) \times 4d + d$
. The Strouhal number (
$St$
) is determined from the dominant frequency of the time series of the spanwise velocity measured in the wake of the array.
Table 1. Long description
A table with 11 rows and 3 columns summarizing the geometric and flow parameters of investigated cases. The columns are labeled Parameter, Symbol, and Value. The table includes parameters such as Stem number, Stem diameter, Stem height, Array width, Solid volume fraction, Frontal area per unit array volume, Stem Reynolds number, Array Reynolds number, Stem flexibility, Stem Strouhal number, and Array Strouhal number. Each row lists the corresponding symbol and value for each parameter. For example, the Stem number is represented by the symbol N with a value of 25, and the Stem diameter is represented by the symbol d with a value of 1.
Instantaneous dynamics of the flexible stem array for different flexibilities (
$K$
), viewed from the side, with (
$a$
)
$K = 1$
, (
$b$
)
$K = 0.32$
, (
$c$
)
$K = 0.1$
and (
$d$
)
$K = 0.032$
. The top-view visualisations highlight the spatial deformation patterns and collective motions of the stems across varying
$K$
. Flow structures are omitted for clarity.
Instantaneous dynamics of the flexible stem array for different flexibilities (
$K$
), viewed from above, with (
$a$
)
$K = 1$
, (
$b$
)
$K = 0.32$
, (
$c$
)
$K = 0.1$
and (
$d$
)
$K = 0.032$
. The top-view visualisations highlight the spatial deformation patterns and collective motions of the stems across varying
$K$
. Flow structures are omitted for clarity.
3.1. Dynamics of the flexible stem array
The dynamics of the stem array is crucial for understanding inter-stem interactions and the resulting flow redistribution within and above the canopy. Figures 2 and 3 present instantaneous snapshots of the array dynamics at different flexibility values (
$K$
) from side and top views, respectively. Distinct dynamic regimes emerge as
$K$
decreases. At
$K = 1$
, the stems behave almost rigidly, exhibiting negligible deformation and oscillations; see figure 2(
$a$
). From the top view (figure 3
$a$
), only minor motion is observed for the downstream stems, confirming the quasi-rigid response of the array. When
$K = 0.32$
, stem deformation becomes more pronounced, as shown in figure 2(
$b$
). As will be shown later, the stem oscillation at this
$K$
is significantly stronger than at other
$K$
values. The top view (figure 3
$b$
) reveals coherent motion across the array, particularly in the downstream rows, where the stems experience stronger fluctuating forces despite the first-row stems directly facing the freestream. At
$K = 0.1$
, the array exhibits patterns similar to those at
$K = 0.32$
, with mild deformation in the first row, and increasingly stronger oscillations downstream, corresponding to intensified vortex activity, as shown in figure 2(
$c$
). The deformation and oscillation amplify towards the array’s trailing edge; see figure 3(
$c$
).
Time evolution of the streamwise tip displacement for stems in different rows, with (
$a$
)
$K = 1$
, (
$b$
)
$K = 0.32$
, (
$c$
)
$K = 0.1$
and (
$d$
)
$K = 0.032$
. For each case, the displacements are ensemble-averaged over all stems within a row.
As
$K$
further decreases, deformation becomes more significant and extends across the entire array. For
$K \leqslant 0.032$
, the first-row stems undergo substantial deformation, indicating that under high flexibility, the leading stems experience significant flow-induced bending, analogous to the lodging state commonly observed in natural aquatic vegetation (Chen & Li Reference Chen, Xu and Ji2025). As shown in figure 2(
$d$
), downstream stems exhibit greater deflection and increasingly complex motions, including higher-order modal dynamics. With weaker structural resistance, stems undergo more pronounced deformation and more irregular motion patterns. Consequently, the flow within and above the array becomes highly heterogeneous. The top-view visualisations (figure 3
$d$
) further confirm extensive deformation across both front and downstream rows, with notably displaced freestream-side stems. These behaviours correspond to intensified vortex interactions and reduced inter-stem shielding, as discussed later.
3.2. Evolution of stem displacements
In this subsection, we examine the temporal evolution of stem displacements in three spatial directions, with particular focus on the streamwise motion of stem tips across the five rows of the array. Figure 4 shows the trajectories of the stem tips over non-dimensional time. The results show that the streamwise displacements increase progressively downstream, and are strongly dependent on the flexibility parameter
$K$
. As
$K$
decreases, the stem tips exhibit larger streamwise excursions from their initial positions, reflecting enhanced structural compliance and hydrodynamic forcing.
With
$K = 1$
, the stems behave almost rigidly, and the maximum non-dimensional streamwise displacement remains below 0.2. As shown in figure 4(
$a$
), the streamwise displacement initially decreases from the first to the second row, and then increases downstream, reaching a secondary maximum near the fourth row. The reduced displacement in the second row results from the shielding effect of the upstream row, which directly faces the incoming flow and induces a substantial momentum deficit in the downstream region. Consequently, the second-row stems experience a diminished drag force and weaker flow excitation. Further downstream, developing vortical structures gradually reenergise the flow, leading to larger displacements. The smaller displacement observed in the last row indicates that for a relatively rigid array, the front row provides the dominant damping on the incoming current.
With
$K = 0.32$
, a similar trend is observed, though the overall deformation amplitude increases. As shown in figure 4(
$b$
), the first-row stems undergo the largest deflection, indicating strong interaction with the freestream and the most effective flow attenuation. The minimum displacement occurs at the fourth row, shifted further downstream compared with the
$K = 1$
case. Notably, oscillatory behaviour emerges in the second and third rows, suggesting enhanced unsteady hydrodynamic coupling among the stems.
When
$K$
is further reduced to 0.1, the streamwise displacement grows markedly, reaching approximately 0.6 in the first row; see figure 4(
$c$
). The pattern of decay and recovery with increasing row number persists, though the oscillations become more pronounced. The lowest displacement is observed in the fourth row, while the streamwise displacements are comparable for the middle three rows. However, with
$K= 0.032$
, a transition occurs: the mid-array stems (R2 and R3) exhibit comparable displacements to those of the first-row stems (see figure 4
$d$
), implying that local drag and flow acceleration within the array govern the deformation more strongly than direct flow exposure.
Across all simulated cases, the results demonstrate that the damping role of successive rows increases once the internal flow develops sufficiently within and above the canopy. This progressive modulation of stem response by flexibility and row position provides key physical insights into momentum absorption, and has important implications for optimising vegetation-based coastal protection and wave attenuation designs.
Time evolution of the vertical tip displacement for stems in different rows, with (
$a$
)
$K = 1$
, (
$b$
)
$K = 0.32$
, (
$c$
)
$K = 0.1$
and (
$d$
)
$K = 0.032$
. The values represent the row-averaged vertical motion of the stem tips.
The temporal evolution of the vertical displacements of stem tips in each row is presented in figure 5. With
$K = 1$
, the vertical displacements are negligible (figure 5
$a$
), consistent with the small streamwise excursions reported earlier, and confirming the minimal deformation of the relatively rigid arrays. As
$K$
decreases, the stems bend progressively downwards, and the vertical displacements become more pronounced (figures 5
$b$
–
$d$
). With
$K = 0.32$
, the downward deflection exceeds that with
$K = 0.1$
, despite the latter being more flexible. This counterintuitive behaviour indicates that the stems exhibiting stronger dynamic responses undergo larger instantaneous deformation. With
$K = 0.032$
, substantial vertical displacements develop, particularly in the middle rows of the array. This pattern mirrors the streamwise deflection shown in figure 4(
$d$
), suggesting that flow acceleration within the array enhances stem bending and reconfiguration. The close correspondence between the vertical and streamwise motions highlights the nonlinear coupling among flexibility, local hydrodynamics and inter-stem interactions, which governs the redistribution of momentum and turbulence within the canopy.
The transverse motion of stem tips is examined next, as vortex-induced vibration (VIV) of flexible stems develops primarily in the cross-flow direction. Figure 6 shows the temporal evolution of transverse tip displacement for arrays with different flexibilities. With
$K = 1$
, the stems remain nearly stationary initially, with oscillations developing at approximately
$t U_{\infty } / d = 15$
. The response is weak in the first column (C1, figure 6
$a$
i), and becomes even smaller towards the array centre (C2 and C3, figures 6
$a$
ii,aiii). In contrast, likely due to reduced downstream confinement and more developed vortex shedding, the stems in the last row exhibit the largest oscillation amplitudes. In the outer columns (C4 and C5), oscillations develop with larger initial amplitudes and slight phase differences (figures 6
$a$
iv,av), implying an early breakdown of flow symmetry.
With
$K = 0.32$
, oscillations develop earlier, and their amplitudes increase markedly. The stems oscillate more synchronously, with larger displacements and well-defined phase relationships, indicating near lock-in between vortex shedding and structural frequency (figures 6
$b$
i–bv). The largest oscillations occur in the outer columns (C1 and C5), whereas the stems in the central region oscillate nearly in phase, highlighting the influence of organised wake structures on array dynamics.
With
$K\leqslant0.1$
, the oscillations become highly irregular, with smaller and more scattered amplitudes (figures 6
$c$
,
$d$
). Motion initiates earlier in the outer columns, whereas the leading stems remain nearly stationary because of strong flow blockage (figures 6
$c$
i,
$c$
v). Downstream stems, exposed to less shielding, exhibit more pronounced oscillations, consistent with less constrained vortex shedding and reduced hydrodynamic coupling between rows. Several downstream stems also display pronounced one-sided deflection (figures 6
$d$
ii,
$d$
iv), suggesting that strong inter-stem interactions and local flow disturbances disrupt the global coherence of the oscillation response.
Time evolution of the transverse tip displacement of the stem array with (
$a$
)
$K = 1$
, (
$b$
)
$K = 0.32$
, (
$c$
)
$K = 0.1$
and (
$d$
)
$K = 0.032$
. Here, S
$i$
denotes the stem numbered
$i$
, as indicated in figure 1.
3.3. Vortex dynamics and inter-stem interactions
Understanding the evolution of vortex dynamics and inter-stem interactions is essential for elucidating the coupling between flow structures and stem responses. Figure 7 shows instantaneous 3-D vortical structures of flexible stem array with different flexibilities, viewed from the top. Distinct flow regimes emerge as
$K$
decreases, highlighting the modulation of wake dynamics by structural flexibility.
With
$K$
= 1, the stems behave nearly rigidly, and the flow field closely resembles that of a rigid array. Vortex shedding initiates at the trailing edge of the first row, with these vortices impinging directly on downstream stems. The separated vortices from subsequent rows are smaller and gradually transform into streamwise vortices. The top view (figure 7
$a$
) reveals a tandem bi-vortex-street pattern (Sumner Reference Sumner2010), with evident lateral interference between adjacent columns. The wake is asymmetric, i.e. narrower on the inner side and broader on the freestream side, reflecting stronger outer vortices and confinement effects (Chang et al. Reference Chang, Constantinescu and Tsai2017). The strongest vortices appear around the third row, after which they dissipate downstream.
With
$K$
= 0.32, the stem response reaches its maximum, and the vortex system becomes considerably more energetic. Alternate vortex shedding occurs from the front-row stems, but due to stem deformation, the vortices shed from the tips are phase-lagged relative to those from the lower portions. The top view (figure 7
$b$
) shows approximately in-phase shedding among the five front stems, followed by strong interactions with newly formed downstream vortices. Beyond the second row, the wake becomes increasingly irregular, and finer vortices fill the inter-column gaps, producing a substantial reduction in streamwise velocity within the array. Above the canopy, the vortices dissipate gradually, without forming large-scale coherent rollers.
With
$K$
= 0.1, the overall vortex intensity weakens slightly. Vortex shedding from the leading stems persists but becomes less organised, with smaller-scale eddies convected downstream. The top view (figure 7
$c$
) shows more irregularity, particularly near the wake of the downstream rows, where direct vortex-stem interactions remain visible. This pattern is reminiscent of the
$K$
= 1 case, indicating a partial recovery of rigid-like behaviour.
With
$K\leqslant 0.032$
, stem deformation is more pronounced, causing out-of-phase shedding and column-dependent differences in wake topology. As shown in figure 7(
$d$
), vortices from the first two rows form narrow wakes, leaving wider gaps for throughflow. Downstream, however, the flow becomes dominated by finer, disordered vortices that occupy the inter-stem spaces. Chaotic motion and asymmetric flow organisation are observed in these large-deformation cases.
The overall patterns of 3-D vortical structures around individual stems are qualitatively similar across all cases, yet discernible differences arise owing to the hydrodynamic interferences among neighbouring stems. Horseshoe vortices consistently develop upstream of the leading-row stems, and to a lesser extent around some downstream stems; see the insets in figure 7(
$a$
). Their strength is, however, noticeably weaker than that observed for a single flexible stem with the same
$K$
(Chang et al. Reference Chang, Constantinescu and Tsai2017; Chang et al. Reference Chen, Ji, Xu, Liu and Campbell2020). This attenuation results from the combined effects of upstream blockage and downstream sheltering within the array; see figures 10 and 13 below.
Instantaneous 3-D vortex dynamics of the flexible stem array, with (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. Vortical structures are visualised using the
$Q$
-criterion (
$Q = 2.5$
; Hunt et al. 1988). Insets of the horseshoe vortex system visualised by the
$x$
-vorticity are provided in (
$a$
). Here,
$\omega _x^{*} \equiv ( {d}/{U_{\infty }}) ({\partial w}/{\partial y} - {\partial v}/{\partial z} )$
.
Instantaneous
$y$
-vorticity fields around and downstream of the flexible stem array at various transverse positions, with (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. The corresponding plane locations are indicated in figure 1. Here,
$\omega _y^{*} \equiv ({d}/{U_{\infty }}) ( {\partial u}/{\partial z} - {\partial w}/{\partial x} )$
.
The influence of stem flexibility on the formation of the horseshoe vortex is also evident. With
$K$
= 1, the stems undergo smaller oscillations, which tend to suppress the near-wall flow separation, and consequently weaken the horseshoe vortex. A similar modulation of the vortex system by stem motion has been reported in the single-stem configuration (Chen & Li Reference Chen, Xu and Ji2025).
In addition to the horseshoe vortex, distinct free-end vortices develop near the stem tips, as revealed by the
$y$
-vorticity distribution in figure 8. These vortices detach into the wake and interact directly with downstream stems. Advected by the upward flow, they convect downstream without approaching the bed. Their subsequent evolution is strongly influenced by the stem dynamics. With
$K$
= 1, where stem deformation and oscillation amplitudes are small, the free-end vortices tend to reattach near the tips of the downstream stems (figure 8
$a$
). As
$K$
decreases, the free-end vortices propagate through the inter-stem gaps and undergo increasingly complex interactions, which are strongly modulated by the transverse stem motions (figures 8
$b$
–
$d$
).
A relatively weak base vortex is also observed behind the trailing stems, corresponding to the lower local velocity induced by the shielding effect. Likewise, the shear layers enveloping the stems are considerably weakened, forming elongated spanwise vortices. The attenuated intensity of these structures reflects the reduced mean velocity within the array. The spanwise vortices play a key role in mediating the exchange of mass and momentum between the canopy and the overlying flow. These spanwise vortices, generated by shear layers at the canopy top, are further modulated by the penetration of outer-layer streamwise vortices. In the wake region, the streamwise vortices emerge from the stretching of spanwise rollers, and subsequently interact with the downstream stems due to the short inter-stem spacing. Furthermore, similar to the single-stem configuration, the hairpin-like vortices are observed in the wake.
The three-dimensionality and intensity of the wake can be quantified through enstrophy, which characterises the spatial distribution of vorticity magnitude (Papaioannou et al. Reference Papaioannou, Yue, Triantafyllou and Karniadakis2006; Thapa et al. Reference Thapa, Zhao, Zhou and Cheng2014). Here, we examine the enstrophy to evaluate how stem flexibility modulates rotational motions within the flow. The enstrophies associated with different components are defined as
\begin{equation} \varepsilon _x=\int _V \omega _x^2\, \text{d} V, \varepsilon _y=\int _V \omega _y^2\, \text{d} V, \varepsilon _{x+y}=\int _V (\omega _x^2+\omega _y^2 )\, \text{d} V, \varepsilon _z=\int _V \omega _z^2\, \text{d} V, \end{equation}
where
$V$
denotes the volume of the computational domain.
Figure 9 shows the variation of the enstrophy components with stem flexibility. As
$K$
decreases, all components exhibit the same non-monotonic trend, increasing to a maximum before decreasing. This behaviour is consistent with the coupling between stem motion and vortical activity, whereby stronger FSI produces larger oscillation amplitudes and more intense vortex dynamics. All enstrophy components reach their maximum at
$K$
= 0.32, corresponding to the strongest stem oscillation. These results indicate that intermediate flexibility promotes the most vigorous 3-D vortex dynamics within the canopy, and consequently the greatest potential for wave-energy dissipation.
Statistics of the vorticity components within the computational domain for flow through flexible stem arrays. The enstrophies associated with different rotational components are compared to illustrate how stem flexibility modulates wake three-dimensionality.
The vertical variation of vortex dynamics provides important insights into the distributions of bed shear stress, turbulence generation and momentum transport within and above the stem array. Figure 10 shows the instantaneous
$z$
-vorticity fields at representative elevations for different flexibility cases. Near the bed (
$z/d = 0.1$
), where the impermeable wall constrains vertical motion, the vortex shedding pattern varies markedly with
$K$
. With
$K$
= 1, distinct vortex shedding occurs from the leading-row stems, but weakens progressively downstream (figure 10
$a$
i), reflecting the reduction in local flow momentum caused by the blockage from the nearly rigid stems. As
$K$
decreases, stem motion at higher elevation exerts a stronger influence on the near-wall flow (figures 10
bi,ci,di). Consequently, vortex shedding behind the stem weakens, whereas vortical activity farther downstream is enhanced, particularly for
$K$
= 0.32, where vigorous oscillation promotes strong local unsteadiness. This modulation of near-wall vortical activity indicates that flexibility modifies both bed shear stress and turbulence generation, the two primary processes governing sediment entrainment (Zhao & Nepf Reference Zhao, Cheng, Wang and Tan2021).
At the mid-canopy elevation (
$z/d = 5.0$
), wall effects become negligible, and the flow is governed by FSI. With
$K$
= 1, two distinct shedding regimes emerge (figure 10
$a$
ii): the freestream-side columns generate strong, detached vortices that convect downstream, whereas vortex shedding within the inner columns is suppressed by wake interference, producing elongated shear layers that reattach to the downstream stems. Such shear-layer reattachment inhibits stem oscillation, and stabilises the local flow. In contrast, with
$K$
= 0.32, pronounced alternating vortices develop within the gaps, and impinge directly on downstream stems (figure 10
$b$
ii), leading to larger oscillation amplitudes and enhanced lateral mixing. The more flexible cases (
$K$
= 0.1 and 0.032) exhibit mixed behaviour: shear-layer reattachment dominates in the upstream rows, whereas vortex–stem interactions become increasingly important downstream (figures 10
$c$
ii,
$d$
ii), indicating a gradual transition in the interaction regime with increasing flexibility.
In the upper canopy region (
$z/d = 9.5$
), free-end effects and stem reconfiguration dominate the vortex dynamics. With
$K$
= 1, vortex shedding weakens progressively towards the stem tips (figure 10
$a$
iii). The inner columns remain dominated by shear-layer reattachment, whereas outer stems continue to shed isolated vortices that induce asymmetric oscillations. With
$K$
= 0.32, coherent vortices persist throughout the upper canopy, and penetrate farther downstream (figure 10
$b$
iii), suggesting enhanced momentum transfer and turbulence generation. In contrast, for the highly flexible cases (
$K$
= 0.1 and 0.032), large stem deflections substantially distort the wake topology. Vortices shed from the upstream stems impinge directly on downstream stems (figures 10
$c$
iii,
$d$
iii), but their strength decreases towards the free ends because flow separation is progressively reduced. Overall, increasing flexibility weakens large-scale vortex coherence while enhancing local turbulent mixing, leading to increasingly anisotropic turbulence and enhanced energy dissipation within the canopy wake.
Instantaneous
$z$
-vorticity fields at different vertical planes for flow through flexible stem arrays, with (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. Here,
$\omega _z^{*} \equiv ({d}/{U_{\infty }}) ( {\partial u}/{\partial x} - {\partial v}/{\partial y} )$
.
3.4. Velocity distribution
The spatial distribution of the streamwise velocity provides key insights into how stem flexibility modifies flow resistance, momentum transport and sediment dynamics within vegetation canopies (Wang et al. Reference Wang, Gualtieri and Huai2022b
). Figure 11 shows the instantaneous streamwise velocity fields at representative transverse planes for cases
$K = 1$
, 0.32, 0.1 and 0.032.
Instantaneous streamwise velocity fields around and downstream of the flexible stem array at various transverse positions, with (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. The corresponding plane locations are indicated in figure 1.
For the relatively rigid case (
$K$
= 1), the stems undergo little deformation, and the flow closely resembles that through a rigid array. At
$y/d = 7$
, the incoming flow is strongly decelerated by the leading-row stems (figure 11
$a$
i). Immediately downstream of each stem, recirculation regions (
$u \lt 0$
) develop due to shear-layer reattachment and periodic vortex shedding. The flow within the array remains at a low velocity, while the flow above the array is accelerated relative to inflow (Nepf Reference Nepf, Ghisalberti, White and Murphy2012; Tschisgale et al. Reference Tschisgale, Löhrer, Meller and Fröhlich2021), reflecting momentum redistribution over the vegetated roughness layer. Downstream of the array, the streamwise velocity gradually recovers, forming alternating low- and high-speed streaks associated with the stem-generated shear layers.
As
$K$
decreases, stem reconfiguration progressively alters the wake topology. With
$K$
= 0.32, pronounced stem deformation and oscillation contract the recirculation regions and enhanced local mixing (figure 11
$b$
i; Wang et al. Reference Wang, Cheng and Ting2022a
,
Reference Wang, Gualtieri and Huaib
). The reduced frontal blockage accelerates the flow through the inter-stem gaps, while the reverse-flow regions behind individual stems are substantially weakened. In contrast, for the more flexible cases (
$K$
= 0.1 and 0.032), the less synchronised stem motion produces irregular flow acceleration between adjacent rows, and a more heterogeneous velocity field within the array (figures 11
$c$
i,
$d$
i). Consequently, the velocity above the canopy is slightly lower than that for
$K$
= 0.32.
At
$y/d = 11$
and 15, the streamwise velocity deficit behind the stems becomes more pronounced for
$K$
= 1, 0.1 and 0.032 (figures 11
$a$
ii,aiii,
$c$
ii,ciii,
$d$
ii,diii), indicating stronger blockage and greater momentum loss within the interior of the array. In contrast, with
$K$
= 0.32, vigorous stem oscillation alleviates the velocity deficit (figures 11
$b$
ii,biii), maintaining a more uniform velocity distribution throughout the canopy.
Within the gaps between adjacent columns, the streamwise velocity locally exceeds the incoming velocity, forming high-momentum corridors that gradually weaken towards the canopy exit (not shown).
Instantaneous streamwise velocity fields around and downstream of the flexible stem array at various streamwise positions, with (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. The corresponding plane locations are indicated in figure 1.
The influence of stem flexibility on the streamwise velocity field is shown in figure 12. Along the first stem row, the streamwise velocity locally exceeds the inflow value, forming narrow jet-like streams (figures 12
$a{-} d$
, row i). Near the stem tips, however, the flow is increasingly governed by stem motion, reflecting the combined effects of flow separation and deformation-induced blockage, which are most pronounced for the most flexible case (
$K$
= 0.032). Downstream, the wakes of the upstream stems strongly modify the local flow, producing alternating high- and low-velocity streaks between successive rows (figures 12
$a{-} d$
, rows ii–v). Their spatial distribution varies markedly with
$K$
, owing to the combined effects of stem deformation and oscillation.
Within the inter-row gaps, the flow is clearly partitioned into low- and high-momentum regions (not shown). Low-momentum zones develop in the stem wake, whereas the high-momentum flow is confined to the inter-stem gaps. Their size and persistence depend strongly on the stem dynamics. For
$K$
= 0.32, pronounced stem oscillations broaden the wake regions and increase their temporal variability, whereas they remain relatively narrow and steady in the nearly rigid case. Farther downstream, the influence of stem motion becomes increasingly pronounced, promoting the wake expansion, and enhancing the spatial heterogeneity of the velocity field. The associated momentum redistribution also contributes to localised high-momentum flow along the outer edge of the array.
Instantaneous streamwise velocity fields around and downstream of the flexible stem array at different vertical elevations, with (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. Streamlines are superimposed as solid lines in each plot.
The distribution of streamwise velocity at different elevations provides crucial insights into the vertical structure of momentum transport, which directly influences sediment entrainment, near-bed shear stress, and the spatial distribution of TKE. Figure 13 shows the instantaneous streamwise velocity fields, superimposed with streamlines, at representative vertical planes for
$K$
= 1, 0.32, 0.1 and 0.032.
Near the impermeable bed (
$z/d = 0.1$
), the streamwise velocity remains small throughout the domain owing to the no-slip boundary condition (figures 13(
$a{-} d$
, row i). The weak shear layers at this elevation correspond to a relatively low local
$Re$
, leading to limited turbulence production. Nevertheless, small-scale velocity fluctuations persist within the array, particularly for
$K$
= 0.32, 0.1 and 0.032, suggesting that stem motion still induces local flow disturbances and weak secondary motions despite strong viscous damping.
At the mid-canopy elevation (
$z/d = 5.0$
), the overall flow pattern is similar across all cases (figures 13
$a$
–
$d$
, row ii). The streamwise velocity reaches its maximum within the inter-column gaps and along the canopy edges, where bypass flow dominates. The internal structure of these high-speed corridors, however, depends strongly on stem flexibility. For
$K$
= 0.32, vigorous stem oscillations periodically constrict the jet-like flow, narrowing the high-speed corridors and enhancing the local shear (figure 13
$b$
ii). By contrast, the remaining cases (
$K$
= 1, 0.1 and 0.032) exhibit weaker velocity modulation (figures 13
$a$
ii,
$c$
ii,
$d$
ii), consistent with their less pronounced stem oscillation.
Near the free ends (
$z/d = 9.5$
), the flow becomes increasingly influenced by stem motion and vertical momentum exchange. For
$K$
= 1, the streamwise velocity within the inter-column gaps increases with height (figure 13
$a$
iii), reflecting flow acceleration above the canopy and the downward entrainment of high-momentum fluid. In contrast, for
$K$
= 0.32, strong oscillations near the stem tips reduce the mean velocity within the array, producing elongated low-velocity regions (figure 13
$b$
iii). For
$K$
= 0.1 and 0.032, the velocity distribution shows intermediate characteristics, i.e. localised acceleration within the gaps, and persistent low-speed wakes behind the stems. Meanwhile, the downstream velocity recovers more rapidly with increasing height, indicating enhanced momentum exchange with the overlying flow, and the gradual transition from the canopy flow to the outer boundary layer.
3.5. Hydrodynamic drag
The influence of stem flexibility on flow damping is manifested through the hydrodynamic drag acting on the stems. The spatial distribution of the streamwise drag force within the array provides valuable insight into the FSI and the overall energy dissipation of vegetated canopies. In this subsection, we present the time-averaged drag force on the stems for the four simulated cases.
Distribution of the mean hydrodynamic drag on each row of stems, with
$K$
= 1, 0.32, 0.1 and 0.032. The drag coefficient for each row is obtained by averaging over five stems.
Figure 14 presents the drag coefficient on the stems in each row for cases
$K$
= 1, 0.32, 0.1 and 0.032. Here, the drag coefficient of each row is defined as
$C_d=({1}/{5}) \sum C_{d, i}=({1}/{5}) \sum F_{d, i} / 0.5 \rho U_{\infty }^2 {d} h$
, where
$F_{d,i}$
denotes the drag force of stem
$i$
. Across all simulated cases, the first row exhibits the largest drag coefficient. This indicates that, regardless of flexibility, the stems at the leading edge experience the strongest resistance to the incoming flow (Koken & Constantinescu Reference Koken and Constantinescu2021; Koken & Constantinescu Reference Koken and Constantinescu2023). This observation is physically reasonable, as the high-velocity incoming flow first impinges on these stems, producing a pronounced deceleration and near-zero velocity at their upstream surfaces.
Downstream of the first row, two distinct trends in the drag coefficient are observed. In case
$K$
= 0.32, the drag coefficient decreases monotonically with increasing row number, suggesting a progressive reduction in flow damping by successive rows. In contrast, for cases
$K$
= 1, 0.1 and 0.032, the drag coefficient attains its minimum at the second row, and subsequently increases gradually, albeit at a slower rate. The gradual increase is closely associated with the development of the gap flow between stems (Koken & Constantinescu Reference Koken and Constantinescu2021).
The FFT spectra of the spanwise velocity and the transverse tip displacement for case
$K$
= 0.32. The velocity is sampled at
$x/d = 35$
and
$y/d = 15$
. The dominant oscillation frequency and the natural frequency are identified from the corresponding FFT spectra.
The variation in drag among the flexible stem arrays indicates that stems with moderate flexibility (case
$K$
= 0.32) exhibit the most effective damping. Figure 15 compares the dominant oscillation frequency, identified from the fast Fourier transform (FFT) of the spanwise velocity, with the natural frequency obtained from the FFT of the transverse tip displacement. The two frequencies coincide for case
$K$
= 0.32, indicating a lock-in state (Rosti et al. Reference Rosti, Banaei, Brandt and Mazzino2018; Rota et al. Reference Rota, Monti, Olivieri and Rosti2024). This finding is in agreement with our previous work (Chen & Li Reference Chen, Xu and Ji2025), which demonstrated that a single flexible stem exhibits its largest oscillation amplitude under lock-in conditions. For both stiffer and more compliant stems (i.e. higher or lower
$K$
), this frequency matching is not achieved, leading to weaker stem responses and reduced flow damping. These results collectively elucidate the roles of flexibility in modulating hydrodynamic drag and energy dissipation within flexible vegetation arrays. Furthermore, the drag coefficients of stems within arrays are markedly lower than those of an isolated flexible element reported by Chen & Li (Reference Chen, Xu and Ji2025). This reduction largely arises from pronounced intra-array shielding, where the mutual sheltering among neighbouring stems attenuates the local velocity and momentum exchange, thereby diminishing the effective drag experienced by individual elements (Etminan et al. Reference Etminan, Lowe and Ghisalberti2017).
3.6. Bed shear stress
The distributions of bed shear stress and TKE are fundamental for predicting sediment transport within vegetated arrays (Yang & Nepf Reference Yang and Nepf2018; Zhao & Nepf Reference Zhao, Cheng, Wang and Tan2021). Gaining a mechanistic understanding of sediment transport within flexible arrays is helpful for advancing coastal protection and ecosystem resilience. Here, the bed shear stress and TKE are defined as
$\tau _b = \rho u_\tau ^2$
and
$k=(\overline {u^{\prime 2}}+\overline {v^{\prime 2}}+\overline {w^{\prime 2}}) / 2$
, respectively, where
$u_{\tau }$
is the local bed friction velocity, and
$u^{\prime }$
,
$v^{\prime }$
and
$w^{\prime }$
are the velocity fluctuations in the streamwise, transverse and vertical directions, respectively. Following Koken & Constantinescu (Reference Koken and Constantinescu2021),
$u_{\tau }$
is calculated using
$u_\tau =(u_m \nu / n_1)^{1 / 2}$
, where
$u_m$
is the velocity magnitude in a plane parallel to the bed situated at a distance
$n_1$
from the wall. The first grid point above the wall is selected to estimate
$u_{\tau }$
. As suggested by Chang et al. (Reference Chang, Constantinescu and Tsai2017) and Cheng et al. (Reference Cheng, Koken and Constantinescu2018), the sediment entrainment capacity of the flow depends on not only the mean near-bed friction velocity (
$\bar {u}_\tau$
) but also the r.m.s. value of the bed friction velocity (
$u_\tau ^{S D}$
), with the latter accounting for the influence of large-scale unsteadiness on the bed shear stress.
Figure 16 shows the mean bed shear stress distribution for cases
$K = 1{-}0.032$
. Overall, both the instantaneous and mean fields exhibit broadly similar spatial patterns across the different flexibility conditions. In all cases, the largest shear stresses occur along the lateral sides of the first-row stems, consistent with previous observations (Chang et al. Reference Chang, Constantinescu and Tsai2017; Chang et al. Reference Chen, Ji, Xu, Liu and Campbell2020). This indicates that, irrespective of whether the stems behave rigidly or flexibly, the strongest bed forcing is concentrated at the upstream edge of the array. Downstream of the first row, however, clear differences emerge that are closely tied to stem flexibility. At intermediate flexibility, where the stem oscillations are most pronounced, the mean bed shear stress is substantially amplified compared with the other cases. This suggests that sediment transport is highly sensitive to stem oscillation amplitude. In addition, a slight increase in bed shear stress is observed near the downstream outlet, attributable to the recovery of the mean streamwise velocity, as is also evident in figures 7, 8 and 10.
Time-averaged bed shear stress (
$\overline {\tau }_b^*$
) at the first grid above the bed, for cases (
$a$
)
$K = 1$
, (
$b$
)
$K$
= 0.32, (
$c$
)
$K$
= 0.1 and (
$d$
)
$K$
= 0.032. The averaging window is not sufficiently long; some instantaneous features remain visible.
Distribution of the r.m.s. bed shear friction velocity at the first grid above the bed for (
$a$
)
$K$
= 1, (
$b$
)
$K$
= 0.32 and (
$c$
)
$K$
= 0.032.
Figure 17 shows the distributions of the r.m.s. bed friction velocity (
$u_{\tau }^{SD}$
) at the first grid above the bed for cases
$K$
= 1, 0.32 and 0.032. In all three cases, elevated values of
$u_{\tau }^{SD}$
occur downstream of the last two stem rows (figures 17
$a$
–
$c$
), corresponding to the formation and downstream evolution of wake vortices. The spatial distribution of
$u_{\tau }^{SD}$
therefore reflects the coupling between stem motion and turbulence generation. As demonstrated in figures 12 and 13, stem oscillation redistributes TKE towards the upper canopy and the overlying flow, reducing the intensity of velocity fluctuations close to the bed. Consequently,
$u_{\tau }^{SD}$
is significantly lower than in the nearly rigid case. This suggests that the enhanced unsteady motions are concentrated in the upper canopy, whereas turbulence penetration towards the bed is suppressed, leading to weaker bed shear fluctuations despite the intensified canopy-scale dynamics (Chen et al. Reference Chen, Ji, Alam, Xu, An, Tong and Zhao2026c
).
4. Discussion: mechanisms that govern the stem dynamics and hydrodynamic forces
In this section, we examine the fundamental flow physics that governs the stem dynamics and the associated hydrodynamic loading within the array. Such analysis provides critical insights into how multiple mechanisms, including shielding, amplification, blockage and oblique-flow effects, as well as delayed flow separation and inter-stem interactions, collectively regulate the magnitude and temporal characteristics of stem motion and the overall hydrodynamic responses.
(
$a$
) Characteristic flow regions surrounding a circular cylinder, highlighting the separated shear layers and the near-wake recirculation bubble. (
$b$
) Flow regimes for two tandem circular cylinders. From left to right: the extended-body regime, the reattachment regime, and the co-shedding regimes.
4.1. Shielding effect
When flow encounters a circular cylinder, the near-wake topology can be divided into several canonical regions: an upstream stagnation zone, a boundary layer developing along the cylinder surface, accelerated flow around the lateral sides, and a separated flow containing a small recirculation bubble; see figure 18(
$a$
) (Zdravkovich Reference Zdravkovich1997). When a second cylinder is placed within the wake of the first, the hydrodynamic loading on the downstream cylinder may be either reduced or enhanced, depending on its position relative to the upstream wake (Sumner Reference Sumner2010). The shielding effect modifies the local dynamic pressure by limiting direct exposure to the freestream. In the low-momentum wake core (region iv), the downstream cylinder experiences reduced stagnation pressure and consequently lower drag, whereas in the high-speed shear layer (region iii), the elevated dynamic pressure results in enhanced drag.
In tandem configurations, the upstream cylinder substantially alters the flow approaching the downstream one, acting as either a turbulence generator or a wake stabiliser (Lee & Basu Reference Lee and Basu1997). Three canonical regimes arise, primarily governed by the spacing ratio
$s/d$
, and influenced by Reynolds number and turbulence intensity (Ljungkrona et al. Reference Ljungkrona, Norberg and Sunden1991; Xu & Zhou Reference Xu and Nepf2004; Sumner Reference Sumner2010); see figure 18(
$b$
). For small spacings (
$s/d \lt1.2{-}1.8$
), the downstream cylinder lies within the vortex-formation region of the upstream cylinder. The upstream shear layers envelop the downstream cylinder without reattachment, and elongated vortices form only behind the downstream body, producing a narrower wake (Lin et al. Reference Lin, Yang and Rockwell2002). The nearly stagnant gap flow results in negative mean drag and weak lift fluctuations, representing an extended-body regime with strong shielding (Alam et al. Reference Alam, Moriya, Takai and Sakamoto2003). For intermediate spacings (
$1.2{-}1.8\lt s/d \lt3.4{-}3.8$
), the upstream shear layers reattach to the downstream cylinder, while vortex shedding remains confined to the downstream body (Zhou & Yiu Reference Zhou and Yiu2006). Although unsteady gap flow may occur (Igarashi Reference Igarashi1981; Kitagawa & Ohta Reference Kitagawa and Ohta2008), the mean drag remains negative with relatively small lift fluctuations, indicating weakened but persistent shielding. For larger spacings (
$s/d \gt3.4{-}3.8$
), both cylinders shed vortices. The downstream cylinder is excited by impinging upstream vortices, enhancing its own shedding and accelerating wake dissipation (Zhou & Yiu Reference Zhou and Yiu2006). In this co-shedding regime, phase-locking causes both cylinders to shed at the same frequency (Zdravkovich Reference Zdravkovich1987; Alam & Zhou Reference Alam and Zhou2007). The mean drag becomes positive and approaches the isolated-cylinder value, accompanied by large force fluctuations due to the elevated dynamic pressure induced by the incoming vortices (Alam et al. Reference Alam, Moriya, Takai and Sakamoto2003).
The shielding effect represents a key mechanism of flow–structure interference that modulates hydrodynamic loading, energy transfer and multiscale flow coherence within arrays of aquatic vegetation. Upstream elements attenuate both the mean velocity and turbulence intensity inside the canopy (figures 10 and 13), causing the flow to transition from the high-speed mixing layer above the canopy to a low-speed, vortex-dominated region characteristic of a von Kármán vortex street within the canopy (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004). In this region, the mean drag on downstream elements is substantially reduced relative to that of an isolated cylinder.
The shielding effect exhibits strong vertical variation, spanning from shear-layer reattachment to direct vortex impingement. Near the bed, where the local velocity is low, the shear layers shed from the front-row stems remain coherent over long streamwise distances. Their persistent coherence strengthens hydrodynamic shielding, suppresses the dynamic pressure on downstream cylinders, and consequently reduces the mean drag well below the isolated-cylinder reference (figures 10
$a$
ii,
$c$
ii). Moving towards the stem tips, the higher local velocity and larger stem oscillations intensify shear-layer roll-up, promoting earlier vortex shedding (figures 10
$b$
ii,biii). In this upper region, downstream cylinders experience strong vortex impingement, which generally elevates their mean drag.
The influence of
$K$
manifests primarily through changes in stem oscillation amplitude and deformation. For moderately flexible stems (
$K$
= 0.32), vigorous oscillations promote vortex impingement over nearly the entire stem length (figure 10
$b$
), thereby altering the form of shielding experienced by downstream elements, from velocity-deficit-driven shielding to shielding characterised by impinging vortices, which can yield elevated mean drag. In contrast, for stiffer and more flexible cases, stem oscillation remains small, and shear-layer reattachment is more commonly observed across the vertical extent, especially around the front-row cylinders (figures 10
$a$
,
$c$
,
$d$
). It should be noted that for pairs of oscillating cylinders, their relative phase may significantly affect the instantaneous and mean hydrodynamic forces (Alam & Zhou Reference Alam and Zhou2007). However, such phase-dependent effects are not considered in the present force-based discussion.
4.2. Delayed separation mechanism
The separation characteristics of bluff bodies play a central role in determining hydrodynamic loading. For a rigid circular cylinder at moderate to high Reynolds numbers, the drag coefficient is governed primarily by the separation angle because the flow detachment determines both the wake width and leeward pressure deficit (Schewe Reference Schewe1983; Zdravkovich Reference Zdravkovich1997; Jiang Reference Jiang2020; Chen & Li Reference Chen, Zhang, Xu, Ji and Ghisalberti2023). Delayed separation narrows the recirculation region, thereby promoting greater momentum entrainment from the freestream, weakening the pressure difference across the body, and ultimately reducing the drag (Achenbach Reference Achenbach1968; Achenbach & Heinecke Reference Achenbach and Heinecke1981; Zdravkovich Reference Zdravkovich1997; Wu et al. Reference Wu, Shen and Chou2004; Jiang Reference Jiang2020). This canonical behaviour provides a useful baseline for understanding flexible stems in turbulent flows.
Unlike rigid cylinders, flexible stems undergo flow-induced deformation that passively reconfigures their geometry. As a stem bends downstream, it becomes increasingly aligned with the local flow (Nepf Reference Nepf, Ghisalberti, White and Murphy2012), weakening the adverse pressure gradient, and delaying separation. The reduced projected area suppresses vortex shedding by prolonging shear-layer development and shifting the effective separation point downstream (Chen & Li Reference Chen, Xu and Ji2025), thereby reducing both mean drag and fluctuating lift (Box et al. Reference Box, Järvelä and Västilä2021; Chen et al. Reference Chen, Jiang and Nepf2026a ).
Within an array, separation delay becomes a collective process governed by the local velocity deficit, wake sheltering, and inter-stem hydrodynamic interactions. Frontal blockage reduces the effective Reynolds number of individual stems, promoting viscous diffusion and delaying separation. Upstream wakes further reduce the incident momentum on downstream stems, alleviating the adverse pressure gradient, increasing the separation angle, and thereby markedly reducing the drag on sheltered stems.
Flexibility induces additional complexity by generating spatially varying reconfiguration along and among stems, leading to spanwise heterogeneous separation behaviour (Chen et al. Reference Chen, Ji, Xu, Dou and Yu2026b ). Regions of strong bending may sustain attached flow, suppressing large-scale vortex formation, and acting as distributed damping that weakens flow-induced forcing across the array. Stem oscillations further modulate the instantaneous curvature and inclination, altering the pressure distribution that governs separation. When the structural natural frequency approaches the vortex-shedding frequency, fluid–structure resonance may occur (Rosti et al. Reference Rosti, Banaei, Brandt and Mazzino2018), promoting spanwise-coherent bending and further delaying separation.
Overall, delayed separation in flexible stem arrays arises from a synergistic interplay between passive reconfiguration, wake sheltering, reduced effective Reynolds number, and fluid–structure resonance. Together, these mechanisms shape the mean and unsteady loading experienced by both natural vegetation and engineered vegetation-mimicking canopies.
4.3. Amplification effect
For an oscillating cylinder, body motion feeds back into the surrounding flow, modifying vortex formation, shedding frequency, and wake topology (Bearman Reference Bearman2011). During VIV, when the oscillation frequency approaches the structural natural frequency, a lock-in regime develops in which vibration amplitude becomes large, and both the mean and fluctuating drag are substantially amplified (Gopalkrishnan Reference Gopalkrishnan1993; Khalak & Williamson Reference Khalak and Williamson1999; Govardhan & Williamson Reference Govardhan and Williamson2000). This amplification results from strong flow–structure coupling that reorganises vortex dynamics and alters the instantaneous pressure distribution around the cylinder.
Cylinder motion continuously modulates the effective angle of attack and separation location. The alternating acceleration and deceleration of the near-wake fluid distort the boundary layer, promoting separation, delaying reattachment, and strengthening recirculation. Elevated instantaneous relative velocity (
$U_R$
) between the flow and the body intensifies vortex development (figure 19), while oscillatory motion injects additional energy into the wake, enhancing turbulence production and Reynolds stresses. The resulting broadened wake increases the pressure deficit and the form drag, producing larger asymmetric pressure fields, and amplifying both the mean and fluctuating forces.
Schematic diagram of the vortex dynamics of (
$a$
) stationary and (
$b$
) oscillating circular cylinders. Here,
$U_R$
is the relative velocity, which is defined as
$U_R=\sqrt {\dot {y}^2+U^2}$
, with
$\dot {y}$
denoting the cylinder transverse velocity, and
$U$
denoting the oncoming flow velocity.
This mechanism underlies the enhanced mean drag observed in aquatic vegetation arrays. For a circular cylinder, the amplification ratio increases approximately with vibration amplitude, with larger amplitudes producing higher mean drag (Bishop & Hassan Reference Bishop and Hassan1964; Sarpkaya Reference Sarpkaya1978, Reference Sarpkaya1995; Khalak & Williamson Reference Khalak and Williamson1999). Flexible aquatic stems undergo flow-induced oscillations that increase local drag and modify the bulk hydraulic resistance of the canopy (Luhar & Nepf Reference Luhar and Nepf2011; Luhar & Nepf Reference Luhar and Nepf2013; Chen et al. Reference Chen, Ji, Alam, Xu, An, Tong and Zhao2026c ). For an individual stem, intermediate flexibility maximises energy dissipation due to resonance-like drag amplification (Chen & Li Reference Chen, Xu and Ji2025), and comparable behaviour is observed in stem arrays, as illustrated in figure 14.
4.4. Blockage effect
The blockage effect refers to the influence of geometric confinement imposed by channel walls, lateral boundaries or neighbouring obstacles on the flow around an immersed body. It becomes significant when a bluff body occupies a substantial portion of a confined flow passage, thereby modifying the local velocity field. First, the reduced flow area accelerates the fluid, particularly within the gap between the body and the sidewall (Maskell Reference Maskell1965). Second, lateral confinement restricts the natural expansion of the wake, causing it to narrow and reattach closer to the body (Zovatto & Pedrizzetti Reference Zovatto and Pedrizzetti2001). Consequently, the vortex shedding frequency typically increases (Zdravkovich Reference Zdravkovich2003). These changes also alter the hydrodynamic loading: drag generally increases because of the higher relative velocity and steeper pressure gradients, while lift is strongly affected by the altered vortex formation (Okamoto & Takeuchi Reference Okamoto and Takeuchi1975). Under strong confinement, the blockage may even delay or suppress boundary-layer separation (Taneda Reference Taneda1964).
Conceptually, blockage can be viewed as an externally imposed confinement effect generated by surrounding boundaries. By modifying the local dynamic pressure and momentum exchange, it influences both the flow structures and the resulting force coefficients. In the present configuration, the sidewall is located relatively close to the vegetation array. It is reasonable to expect that the confinement accelerates the local flow adjacent to the array, and enhances vortex dynamics in the outer-array region (i.e. between the array edge and the sidewall), as discussed in
$\S$
3.3. These effects produce slightly stronger stem motions and modestly higher hydrodynamic loading.
4.5. Oblique flow effect
The inclination of a circular cylinder relative to the oncoming flow profoundly alters its vortex dynamics, wake topology and hydrodynamic loading. As illustrated in figure 20(
$a$
), the incoming flow can be decomposed into a component normal to the cylinder axis and an axial component aligned with the cylinder axis. The normal component drives the formation of spanwise vortices, whereas the axial component disrupts their coherence (Ramberg Reference Ramberg1983; Kawamura & Hayashi Reference Kawamura and Hayashi1995; Lucor & Karniadakis Reference Lucor and Karniadakis2003; Zhou et al. Reference Zhou, Razali, Zhou, Chua and Cheng2009). Consequently, vortex shedding becomes fully 3-D and oblique, with spanwise phase differences producing helical or skewed vortices (Shirakashi et al. Reference Shirakashi, Hasegawa and Wakiya1986; Wang et al. Reference Wang, He, Dey and Fang2019). The enhanced streamwise vorticity generated by the axial flow promotes vortex stretching, reorientation, and more frequent vortex dislocations arising from non-uniform phase distributions (Lucor & Karniadakis Reference Lucor and Karniadakis2003; Zhao Reference Zhao2015). These processes weaken the spanwise coherence, and diminish the regularity of the von Kármán vortex street. Because the effective cross-flow velocity decreases with inclination, the shedding frequency generally decreases (Van Atta Reference Van Atta1968; Shirakashi et al. Reference Shirakashi, Hasegawa and Wakiya1986), accompanied by reduced mean drag and less periodic force fluctuations.
Schematics of the flow past an inclined (
$a$
) rigid circular cylinder and (
$b$
) flexible vegetation element. The incoming flow approaches from left to right. The inclination angle
$\alpha$
denotes the angle between the freestream direction and the cylinder (or stem) axis.
When the inclined cylinder undergoes VIV, the weaker vortices and reduced coherence attenuate the fluctuating lift, and disrupt spanwise synchronisation, leading to desynchronised or partially locked-in responses. For a flexible vegetation element, however, the dynamics becomes more complex. As depicted in figure 20(
$b$
), the effective inclination angle varies along the stem due to curvature-dependent deformation, producing vertically non-uniform FSIs. The effects of deflection and oscillation for an individual flexible stem have been discussed in Chen & Li (Reference Chen, Xu and Ji2025). Towards the stem tip, increasing deformation produces a more streamlined configuration, while the oscillation amplitude also increases. Consequently, vortex shedding is governed by the combined effects of deformation and oscillation. Under strong reconfiguration, deformation dominates, and vortices are shed primarily from the tip, approximately aligned with the mean flow. When the stem behaves more rigidly, oscillation dominates, and vortex shedding intensifies towards the tip due to the larger oscillation amplitude. More generally, a flexible inclined cylinder develops highly 3-D, weakly coherent vortex structures because of the combined effects of axial flow and structural deformation, resulting in lower effective cross-flow velocity, reduced drag and lift, and suppressed flow-induced vibration (Chen et al. Reference Chen, Jiang and Nepf2026a
,
Reference Chen, Ji, Xu, Dou and Yub
). These inclination-dependent mechanisms are essential for predicting hydrodynamic loading, wake characteristics and structural responses in engineering and environmental applications.
5. Conclusions
In this study, we investigated the dynamics of flow through a square array of flexible stems with varying flexibility using direct numerical simulations. We analysed the stem response, 3-D vortex dynamics and interactions, and flow distributions within and above the stem arrays. Moreover, we discussed the hydrodynamic drag on the stems, and the distributions of wall shear stress and turbulent kinetic energy. The principal findings of this study can be summarised as follows.
The dynamics of stem arrays is strongly influenced by stem flexibility. Three distinct response regimes are identified within the simulated parameter space. At high
$K$
, the stems behave similarly to rigid cylinders, exhibiting only limited deformation and oscillations. At intermediate
$K$
, the structural natural frequency matches the oscillation frequency, producing a lock-in regime with the strongest stem response. At low
$K$
, stem oscillation becomes moderate, while deformation increases markedly. The array dynamics is more irregular, with higher-mode responses emerging. The evolution of stem displacement further suggests that stem deformation generally increases as
$K$
decreases.
The vortex dynamics within and around the stem arrays reveals several interaction mechanisms that explain these response regimes. The intensity of vortical activity reaches its maximum for case
$K$
= 0.32, corresponding to the largest oscillation amplitude. Both stem oscillation and deformation influence hydrodynamic interactions, which can be broadly classified into two categories: shear-layer reattachment, and vortex–vortex interaction. Shear-layer reattachment promotes pressure balance around the stem circumference and suppresses oscillation, while vortex–vortex interactions enhance vortex shedding, intensify pressure asymmetry, and consequently amplify the stem response.
Flow distributions from multiple perspectives further elucidate the hydrodynamics within the stem arrays. The flow accelerates through the gap between adjacent stem columns, forming jet-like structures, whereas substantially lower velocities develop in the wake of the stems because of shielding. Both stem oscillation and deformation strongly influence the flow structure and vortex dynamics. Stem oscillation increases the flow velocity in the upper region of the array, while reducing the near-bed flow momentum. Stem deformation modifies the spacing between adjacent stems, leading to a vertical transition of the dominant interaction mechanism from vortex–vortex interaction to shear-layer reattachment.
The energy-dissipation and sediment-transport capacities of the stem arrays are evaluated using the hydrodynamic drag together with the spatial distributions of the mean bed shear stress (
$\overline {\tau }_b^*$
) and root mean square bed friction velocity (
$u_{\tau }^{SD}$
). The array with intermediate
$K$
exhibits the greatest energy dissipation, surpassing even that of the rigid configuration. The mean bed shear stress reaches its maximum near the upstream edge of the array, and remains appreciable towards the downstream end. Although stem oscillation enhances vortex dynamics, it weakens the near-bed flow, resulting in a lower
$u_{\tau }^{SD}$
than that of the relatively rigid cases.
A unified physical framework emerges in which shielding, delayed separation, drag amplification, blockage and oblique-flow effects collectively govern the dynamics of flexible stem arrays and the resulting hydrodynamic loading. By linking flow topology, wake interference and deformation-induced reconfiguration, this framework explains how stem motion and force generation arise from the interplay between velocity-deficit shielding, resonance-like drag amplification, confinement-modified wake dynamics, and inclination-induced 3-D vortex shedding. These coupled mechanisms collectively contribute to determining the mean and unsteady hydrodynamic forces within vegetated canopies, thereby regulating energy dissipation, momentum transport and structural resistance in aquatic environments.
This study provides physical insight into the regime-dependent flow–structure interactions governing hydrodynamic forces, energy dissipation and transport processes within flexible stem arrays. These results further suggest that arrays with intermediate flexibility may provide an effective strategy for enhancing energy dissipation in coastal protection applications.
The present study considers only a regularly spaced, aligned stem array, whereas natural vegetation canopies generally exhibit irregular and spatially heterogeneous distributions. Alternative stem arrangements may alter wake interactions and the formation of high-velocity corridors, thereby affecting vortex lock-in and inducing subtle changes in the near-bed flow.
Acknowledgements
This work was supported with computational resources provided by the National Supercomputing Centre, Singapore (https://www.nscc.sg), under project IDs 11003759 and 11004188.
Funding
This project is funded by the Research, Innovation and Enterprise 2025 Coastal Protection and Flood Management Research Programme. The authors are grateful to the National Research Foundation of Singapore, PUB (Singapore’s National Water Agency) for supporting this work done in the Coastal Protection and Flood Resilience Institute (CFI) Singapore under the project grant H1-P3.
Declaration of interests
The authors report no conflict of interest.
a
b
4d
h/d=10
D
D=(N−1)×4d+d
St
K
a
K=1
b
K=0.32
c
K=0.1
d
K=0.032
K
K
a
K=1
b
K=0.32
c
K=0.1
d
K=0.032
K
a
K=1
b
K=0.32
c
K=0.1
d
K=0.032
a
K=1
b
K=0.32
c
K=0.1
d
K=0.032
a
K=1
b
K=0.32
c
K=0.1
d
K=0.032
i
i
a
K
b
K
c
K
d
K
Q
Q=2.5
x
a
ωx∗≡(d/U∞)(∂w/∂y−∂v/∂z)
y
a
K
b
K
c
K
d
K
ωy∗≡(d/U∞)(∂u/∂z−∂w/∂x)
z
a
K
b
K
c
K
d
K
ωz∗≡(d/U∞)(∂u/∂x−∂v/∂y)
a
K
b
K
c
K
d
K
a
K
b
K
c
K
d
K
a
K
b
K
c
K
d
K
K
K
x/d=35
y/d=15
τ¯b∗
a
K=1
b
K
c
K
d
K
a
K
b
K
c
K
a
b
a
b
UR
UR=y˙2+U2
y˙
U
a
b
α
Open access
