1. Introduction
Vegetation-generated turbulence influences many processes in river ecosystems. It enhances nutrient uptake rate and promotes plant growth by thinning the diffusive sublayer on individual leaves, as demonstrated by Cornacchia et al. (Reference Cornacchia, Licci, Nepf, Folkard, van der Wal, van de Koppel and Bouma2019b), who found a strong correlation between ammonium uptake rates and total kinetic energy at the channel scale. However, some studies have noted a decrease in plant growth rate at higher turbulence levels (Asaeda & Rashid Reference Asaeda and Rashid2017). Vegetated turbulence also impacts the suitability for larval habitat (e.g. Prada et al. Reference Prada, George, Stahlschmidt, Jackson, Chapman and Tinoco2021; Tinoco et al. Reference Tinoco, Prada, George, Stahlschmidt, Jackson and Chapman2022). With respect to fluvial morphodynamics, vegetation-generated turbulence can promote sediment resuspension and increase sediment transport rates, relative to unvegetated channels with the same velocity (e.g. Tinoco & Coco Reference Tinoco and Coco2016; Yang & Nepf Reference Yang and Nepf2018; Shan et al. Reference Shan, Zhao, Liu and Nepf2020; Zhao & Nepf Reference Zhao and Nepf2021). This has led to the development of predictive models for bedload transport based on channel-averaged turbulent kinetic energy (TKE), rather than mean velocity (e.g. Yang & Nepf Reference Yang and Nepf2019; Shan et al. Reference Shan, Zhao, Liu and Nepf2020). Additionally, vegetation impacts turbulent diffusion throughout the water column, altering the vertical profile of suspended sediment concentration (Xu et al. Reference Xu and Nepf2020) and enhancing gas transfer rates across both the air–water and sediment–water interfaces (Tseng & Tinoco Reference Tseng and Tinoco2020, Reference Tseng and Tinoco2022). Because numerous processes in natural streams are impacted by vegetation-generated turbulence and various models for them have been developed based on channel-averaged TKE, it is important to be able to predict the levels of turbulence in channels and to understand how turbulence production is impacted by changes in plant area density and spatial configuration. These are the goals of the present study.
Figure 1. (a) Schematic of 3-D flow structure around a canopy of 3-D submerged patches of model vegetation. Black arrows illustrate the time-mean velocity profiles. Red and blue curves illustrate turbulent wake vortices created in the horizontal (red) and vertical (blue) plane, respectively. Vertical distribution of time- averaged and horizontally averaged (b) streamwise velocity and (c) RS within the canopies shown in side view. Here
$H$
is the water depth,
$h$
is the deflected canopy height,
$\tilde{h}$
is the location of maximum RS,
$L_{p}$
is the patch length and
$\theta$
is the inclination angle of the patch;
$U_{A}$
and
$U_{v}$
denote the average velocity above and within the vegetation, respectively, and
$\tilde{U}_{A}$
and
$\tilde{U}_{B}$
are the mean velocity above and below
$\tilde{h}$
, respectively;
$\delta _{e}$
is the vortex penetration length scale and
$\delta _{s}$
is the distance from the canopy top to
$\tilde{h}$
(
$\delta _{s}=h-\tilde{h}$
).
Natural streams are often inhabited by submerged vegetation patches with finite length and width, hereafter called three-dimensional (3-D) submerged patches (Sukhodolova & Sukhodolov Reference Sukhodolova and Sukhodolov2012; Schoelynck et al. Reference Schoelynck, Creëlle, Buis, De Mulder, Emsens, Hein and Folkard2018; Yan et al. Reference Yan, Zhu, Zhou, Chu, Sui, Cui and van der Heide2021; Calvani, Carbonari & Solari Reference Calvani, Carbonari and Solari2022). When flow encounters an individual 3-D submerged patch, a turbulent wake is produced with shear layers formed around and at the top of the patch (figure 1 a). In a natural stream with numerous 3-D patches, turbulent wakes of individual patches interact with surrounding ones, leading to the development of a complex canopy flow (Biggs et al. Reference Biggs, Nikora, Gibbins, Cameron, Papadopoulos, Stewart and Hicks2019; Cornacchia et al. Reference Cornacchia, Folkard, Davies, Grabowski, van de Koppel, van der Wal and Bouma2019a ). Three different flow regimes have been observed, depending on the size of the gap between patches (Folkard et al. Reference Folkard2011; Mayaud, Wiggs & Bailey Reference Mayaud, Wiggs and Bailey2016): skimming flow, wake interference flow and isolated roughness flow. Skimming flow occurs when patches are spaced close enough that distinct wakes cannot develop, and the patchy canopy acts hydrodynamically as a single continuous canopy. In this regime, most of the flow skims over the top of the canopy. In the second regime, wake interference flow, the patch spacing is larger, allowing individual wakes to develop, but wakes between patches partially interact. In this regime, flow accelerates both between patches and over the top of the patchy canopy. Finally, isolated roughness flow occurs when patches are distributed widely enough that individual wakes do not interact with each other.
Many previous studies have investigated the turbulent flow characteristics in the gap between submerged patches across a wide range of spatial densities (e.g. Folkard Reference Folkard2005, Reference Folkard2011; Maltese et al. Reference Maltese, Cox, Folkard, Ciraolo, La Loggia and Lombardo2007; Chung et al. Reference Chung, Mandel, Zarama and Koseff2021). These studies have primarily focused on the mean and turbulent properties within the gap (e.g. Folkard Reference Folkard2005, Reference Folkard2011; Chung et al. Reference Chung, Mandel, Zarama and Koseff2021) or on the turbulent coherent structures at the canopy–gap interface (Maltese et al. Reference Maltese, Cox, Folkard, Ciraolo, La Loggia and Lombardo2007). While these studies defined important characteristics of turbulent flow within the gap, they only considered channel-spanning patches, and were thus limited to a two-dimensional (2-D) flow field. Some studies have considered 3-D patches, but they focused on evaluating hydraulic resistance using the Manning’s coefficient in different spatial patch configurations (e.g. Green Reference Green2005; Savio et al. Reference Savio, Vettori, Biggs, Zampiron, Cameron, Stewart and Nikora2023) or analysing the patch growth pattern based on the drag measured around the patch (e.g. Cornacchia et al. Reference Cornacchia, Folkard, Davies, Grabowski, van de Koppel, van der Wal and Bouma2019a ), with little consideration of turbulent flow characteristics. Therefore, there is a need for deeper analysis of the 3-D velocity field and turbulent structures associated with 3-D submerged patches.
The aim of this study was to investigate how the spatial configuration of 3-D submerged patches impacts flow characteristics and turbulence production, and to develop a predictive model for the channel-averaged TKE. Laboratory experiments were conducted using patches of model plants. Both longitudinal and lateral spacing between patches was varied to create a range of patch configurations that included both skimming flow and wake interference flow. Detailed velocity measurements were used to analyse the sources of turbulence production within a submerged canopy: (i) wake production and (ii) canopy shear production. A physically based predictive model for the channel-averaged TKE was developed by extending previous TKE models developed in homogenously distributed canopies (Tanino & Nepf Reference Tanino and Nepf2008; Xu & Nepf Reference Xu and Nepf2020; Zhang et al. Reference Zhang, Lei, Huai and Nepf2020), which was verified by comparison with measurements.
2. Predictive model for channel-averaged TKE
This study developed a prediction for channel-averaged TKE in a channel with a submerged canopy of 3-D vegetation patches (figure 1). Wilson & Shaw (Reference Wilson and Shaw1977) introduced the TKE equation in canopy flow based on temporal averaging and spatial averaging schemes, removing small-scale turbulence and spatial heterogeneity at the scale of individual plant elements. As described in detail in Raupach & Shaw (Reference Raupach and Shaw1982), under steady and fully developed flow conditions, the double-averaged TKE budget is
\begin{equation}{{\partial {k_t}} \over {\partial t}} = {P_w} + {P_s} + T - \varepsilon = 0,\end{equation}
with
\begin{equation} P_{w}=-\left\langle \overline{u^{\prime}_{i}u^{\prime}_{j}}^{\prime\prime}\frac{\partial \overline{u}^{\prime\prime}_{i}}{\partial x_{j}}\right\rangle ,\end{equation}
\begin{equation} P_{s}=-\big\langle \overline{u^{\prime}w^{\prime}}\big\rangle \frac{\partial \langle \overline{u}\rangle }{\partial z}, \end{equation}
\begin{equation} T=\frac{\partial }{\partial z}\left(-\frac{\big\langle \overline{u^{\prime}_{i}u^{\prime}_{i}w^{\prime}}\big\rangle }{2}-\frac{\big\langle \overline{p^{\prime}w^{\prime}}\big\rangle }{\rho }+\nu \left\langle \overline{u^{\prime}_{i}\left(\frac{\partial u^{\prime}_{i}}{\partial z}+\frac{\partial w^{\prime}}{\partial x_{i}}\right)}\right\rangle -\frac{\big\langle \overline{u^{\prime}_{i}u^{\prime}_{i}}^{\prime\prime}\overline{w}^{\prime\prime}\big\rangle }{2}\right). \end{equation}
Here,
$k_{t}$
is the TKE per fluid mass,
$P_{w}$
,
$P_{s}$
and
$\varepsilon$
are the wake production, shear production and turbulent dissipation rate, respectively, and
$T$
is the transport term, containing turbulent, pressure, viscous and dispersive transport terms, shown in order in (2.4). Subscripts i and j are equal to 1, 2 and 3, corresponding to the streamwise (x), lateral (y) and vertical (z) directions, respectively. The overbar and angle brackets indicate a temporal and spatial average, and a single and double prime signifies temporal and spatial fluctuations, respectively. The spatial average (bracket) is performed over thin horizontal slabs that are sufficiently wide to remove the heterogeneity at the canopy-element scale but thin enough to preserve vertical gradients, and only includes the fluid domain, excluding the volume occupied by canopy elements (Raupach & Shaw Reference Raupach and Shaw1982).
Averaging (2.1) over the flow depth
$H$
, the transport terms drop out. While the dispersive transport term within the canopy contributes to turbulence production, it accounts for only 2.5 % of the total channel-averaged turbulence in the sparsest configuration with the largest spatial variations in flow characteristics. Therefore, we exclude this term from our channel-averaged TKE model. Thus, for steady flow conditions, the channel-averaged turbulent dissipation rate balances with the summation of wake production and shear production:
\begin{equation} \left\langle \varepsilon \right\rangle _{z}=\left\langle P_{w}\right\rangle _{z}+\left\langle P_{s}\right\rangle _{z}. \end{equation}
Here, 
denotes the average over both width (bracket) and depth (subscript ‘z’). While the shear production
$\langle P_{s}\rangle _{z}$
is typically neglected for emergent canopies with negligible vertical shear (e.g. Nepf & Vivoni Reference Nepf and Vivoni2000; Xu & Nepf Reference Xu and Nepf2020), this term must be considered for submerged canopies because strong turbulence is produced at the top of the canopy (e.g. Poggi, Katul & Albertson Reference Poggi, Katul and Albertson2004; King, Tinoco & Cowen Reference King, Tinoco and Cowen2012).
Turbulent wakes generated by individual plants or patches of plants convert mean flow kinetic energy into TKE. If we assume that there is no energy loss in the conversion process, the wake turbulence production is equal to the mean flow energy dissipation by canopy form drag. Based on this, Raupach & Shaw (Reference Raupach and Shaw1982) represented the cylinder-wake production as the rate of work done by form drag,
$P_{w}=-\langle \overline{u^{\prime}_{i}u^{\prime}_{j}}^{\prime\prime}({\partial \overline{u}^{\prime\prime}_{i}}/{\partial x_{j}})\rangle \approx 0.5C_{D}^{form}a\langle \overline{u}\rangle ^{3}$
, from which
\begin{equation} \left\langle P_{w}\right\rangle _{z}=\frac{1}{H}\int _{0}^{H}\frac{1}{2}C_{D}a\langle \overline{u}\rangle ^{3}{\rm d}z\approx \frac{C_{D}aL_{c}U_{v}^{3}}{2H}, \end{equation}
in which
$C_{D}^{form}$
is the form drag coefficient,
$a$
is the vegetation frontal area per canopy volume and
$H$
is the water depth. The form drag coefficient was assumed to be equal to the total drag coefficient (
$C_{D}^{form}\approx C_{D}$
) because the viscous drag generated by the wetted surface area of individual stems/leaves makes a negligible contribution to wake turbulence production. Specifically, for the patches with high internal solid volume fraction, the drag is dominated by the shear layers that form around the patch, rather than internal viscous and form drag on individual elements within the patch (Nicolle & Eames Reference Nicolle and Eames2011). Here
$L_{c}$
is a length scale over which the plant presents the frontal area toward the flow and, thus, generates drag. Note that natural vegetation is flexible and can be deflected in response to the flow, called reconfiguration, which alters the height of the patch and the orientation of the leaves along the stems, which can impact the vegetation drag (Sand-Jensen Reference Sand-Jensen2003; Siniscalchi & Nikora Reference Siniscalchi and Nikora2013; Vogel Reference Vogel2020). For plants consisting of individual stems (e.g. reeds and sedges), the frontal area is only contributed by the vertical component of the deflected stem or leaf, such that
$L_{c}$
is the deflected canopy height
$h$
(e.g. King et al. Reference King, Tinoco and Cowen2012; Lei & Nepf Reference Lei and Nepf2021). For plants with leaves distributed along a stem (e.g. Rotala indica, figure 1; Potamogeton natans, Vettori et al. Reference Vettori, Niewerth, Aberle and Rice2021; Ranunculus penicillatus, Biggs et al. Reference Biggs, Nikora, Gibbins, Cameron, Papadopoulos, Stewart and Hicks2019), the leaves along the deflected part of the stem contribute to the frontal area and drag (see also Zhang & Nepf Reference Zhang and Nepf2020), so that
$L_{c}$
is larger than
$h$
. For this morphology, the frontal area length scale
$L_{c}$
can be estimated geometrically from the inclination angle
$\theta$
describing the degree of patch reconfiguration (see figure 1
c). Specifically, by simplifying the geometry of the deflected patch as a right triangle with a base of
$L_{p}$
and a perpendicular height of
$h$
(figure 1
c), the drag acting along the deflected part (
$F_{D,L}$
) can be approximated as
$F_{D,L}=F_{D,h}\cos \theta$
(yellow arrow in figure 1
c), in which
$F_{D,h}$
is the drag acting on the perpendicular part. Adjusting the direction of
$F_{D,L}$
to oppose the flow by multiplying by
$\cos \theta$
, the total drag acting on the deflected patches is
$F_{D,\textit{total}}=F_{D,h}(1+\cos ^{2}\theta )$
, leading to
$L_{c}=h(1+\cos ^{2}\theta )=h(1+ {L_{p}^{2}/}{h^{2}+L_{p}^{2}})$
, with
$L_{p}$
defining the patch length in the horizontal plane. With this definition, the frontal area per canopy volume (
$a$
) was computed as
$a=n_{\textit{total}}W_{p}L_{c}/L_{test}WH$
, in which
$n_{\textit{total}}$
is the total number of patches, each of width
$W_{p}$
, that are distributed over a channel of width
$W$
and a test section of length
$L_{test}$
. Here
$U_{v}$
is the mean velocity in the vegetated layer (
$0\lt z\lt h$
), which can be evaluated from the momentum balance (e.g. Nepf Reference Nepf2012; Luhar & Nepf Reference Luhar and Nepf2013):
\begin{equation} \rho gSH-\tau _{b}-\tfrac{1}{2}\rho C_{D}aL_{c}U_{v}^{2}=0. \end{equation}
Here, S is the water surface slope, g is the gravitational acceleration,
$\rho$
is the density of water,
$\tau _{b}$
is the bed shear stress.
Assuming canopy drag is large compared with the bed friction (
$\tau _{b}\ll 0.5\rho C_{D}aL_{c}U_{v}^{2}$
), (2.7) predicts that
\begin{equation} U_{v}=\sqrt{\frac{2gSH}{C_{D}aL_{c}}}, \end{equation}
and combining (2.6) and (2.8), the channel-averaged wake production is
\begin{equation} \left\langle P_{w}\right\rangle _{z}=gS\sqrt{\frac{2gSH}{C_{D}aL_{c}}}. \end{equation}
The assumption of
$\tau _{b}\ll 0.5\rho C_{D}aL_{c}U_{v}^{2}$
was justified based on our measurement, where the range of bed shear stress (
$\tau _{b}=\rho C_{f}U_{v}^{2}$
) was O(0.0001–0.001), which was negligible compared with the canopy drag with O(0.01–0.1). Here, the bed-friction coefficient was calculated as
$C_{f}=-\overline{u^{\prime}w^{\prime}}_{nb}/U_{0}^{2}=0.003$
, in which
$\overline{u^{\prime}w^{\prime}}_{nb}$
is the near-bed Reynolds stress (RS) measured over bare bed and
$U_{0}$
is the channel-averaged streamwise velocity.
Next, consider the shear production term. For a uniform array of submerged stems, Zhang et al. (Reference Zhang, Lei, Huai and Nepf2020) developed a prediction for canopy-average shear production that assumed the peak RS occurred at the top of the canopy, which is typically observed for canopies consisting of a uniform array of cylinders or plant models (e.g. Poggi et al. Reference Poggi, Katul and Albertson2004; Chung et al. Reference Chung, Mandel, Zarama and Koseff2021). However, measurements made in the present study revealed that in a patchy canopy the location of the maximum RS was influenced by the distance between patches within the canopy. When the patches were closely spaced, the maximum RS was at the top of the canopy, but the position dropped below the top of the canopy as the spacing between patches increased. We denote the location of the maximum RS as
$z=\tilde{h}$
(figure 1
c), and the RS at
$z=\tilde{h}$
can be described as in Chen, Jiang & Nepf (Reference Chen, Jiang and Nepf2013),
\begin{equation} \big\langle -\overline{u^{\prime}w^{\prime}}\big\rangle _{max}=\big\langle -\overline{u^{\prime}w^{\prime}}\big\rangle _{z=\tilde{h}}=u_{*,z=\tilde{h}}^{2}=C\big(\tilde{U}_{A}-\tilde{U}_{B}\big)^{2}, \end{equation}
with the empirical coefficient
$C$
, based on scaling arguments in Gioia & Bombardelli (Reference Gioia and Bombardelli2002),
\begin{equation} C=K_{c}\left(\frac{\delta _{e}}{H}\right)^{1/3}, \end{equation}
in which
$u_{*,z=\tilde{h}}$
is the shear velocity at
$z=\tilde{h}$
,
$\tilde{U}_{A}$
and
$\tilde{U}_{B}$
are the mean velocity above and below
$\tilde{h}$
, respectively (figure 1
b), and
$K_{c}$
is an empirical factor depending on the vegetation geometric characteristics, which will be determined in the results section. The length scale
$\delta _{e}$
characterizes the penetration of turbulent momentum flux into the canopy, which is defined as the distance between the point of the maximum RS and the point where the RS has decayed to 10 % of its maximum value (Nepf & Vivoni Reference Nepf and Vivoni2000). For homogeneously distributed canopies with
$C_{D}ah\gt 0.1$
, the location of the maximum RS is at the top of the canopy (Nepf et al. Reference Nepf, Ghisalberti, White and Murphy2007). However, for a patchy canopy, the peak in RS is shifted down from the canopy top (
$z=h$
) by a distance
$\delta _{s}$
(i.e.
$\delta _{s}=h-\tilde{h}$
, figure 1
c). This study experimentally determined the variables
$\delta _{s}$
,
$\delta _{e}$
and
$K_{c}$
, which will be described in the result section.
Ghisalberti (Reference Ghisalberti2009) demonstrated that the difference between the velocity at the top of a canopy,
$U_{h}$
, and the velocity within the canopy,
$U_{v}$
(given by (2.8)), is related to the shear velocity at the top of the canopy. Specifically,
$U_{h}-U_{v}\approx 2.6u_{*,z=h}$
. Assuming this holds for the peak RS, even as it is shifted downward to
$z=\tilde{h}$
for a patchy canopy,
\begin{equation} \tilde{U}_{B}\approx U_{\tilde{h}}-2.6u_{*,z=\tilde{h}}=U_{\tilde{h}}-2.6\sqrt{gS\big(H-\tilde{h}\big)}, \end{equation}
in which
$U_{\tilde{h}}$
is time-mean velocity at
$z=\tilde{h}$
.
The velocity shear at
$z=\tilde{h}$
can be predicted as (e.g. Zhang et al. Reference Zhang, Lei, Huai and Nepf2020),
\begin{equation} \frac{\partial \langle \overline{u}\rangle }{\partial z}_{z=\tilde{h}}\approx \frac{U_{\tilde{h}}-\tilde{U}_{B}}{\delta _{e}}=\frac{2.6u_{*,z=\tilde{h}}}{\delta _{e}}=\frac{2.6}{\delta _{e}}C^{1/2}\big(\tilde{U}_{A}-\tilde{U}_{B}\big). \end{equation}
Substituting (2.10) and (2.13) into (2.3), the shear production at
$z=\tilde{h}$
is
\begin{equation} P_{s,z=\tilde{h}}=-\overline{u^{\prime}w^{\prime}}_{z=\tilde{h}}\frac{\partial \langle \overline{u}\rangle }{\partial z}_{z=\tilde{h}}=\frac{2.6}{\sqrt{H\delta _{e}}}K_{c}^{\frac{3}{2}}\big(\tilde{U}_{A}-\tilde{U}_{B}\big)^{3}, \end{equation}
with
\begin{equation} \tilde{U}_{A}=\frac{U_{0}H-\tilde{U}_{B}\left(h-\delta _{s}\right)}{H-\left(h-\delta _{s}\right)}. \end{equation}
Within the canopy, shear production decreases from a peak at
$z=\tilde{h}$
to zero at the penetration depth
$z=\tilde{h}-\delta _{e}$
. Similarly, shear production above the canopy approaches to zero at the water surface
$z=H$
. For simplicity, we can assume a linear decrease of shear production within and above the canopy, such that the channel-averaged shear production
$\langle P_{s}\rangle _{z}$
becomes
\begin{equation} \left\langle P_{s}\right\rangle _{z}=-\frac{1}{H}\int _{0}^{H}\big\langle \overline{u^{\prime}w^{\prime}}\big\rangle \frac{\partial \langle \overline{u}\rangle }{\partial z}\mathrm{d}z=\frac{1.3}{H\sqrt{H\delta _{e}}}K_{c}^{\frac{3}{2}}\big(\tilde{U}_{A}-\tilde{U}_{B}\big)^{3}\left(H-h+\delta _{s}+\delta _{e}\right). \end{equation}
The turbulent dissipation rate scales with a ratio of TKE to characteristic eddy scale,
$L_{e}$
(Tennekes & Lumley Reference Tennekes and Lumley1972):
\begin{equation} \varepsilon \sim \frac{k_{t}^{3/2}}{L_{e}}. \end{equation}
Within an emergent canopy,
$L_{e}$
is equal to the canopy-element scale, e.g. stem diameter or leaf width (e.g. Tanino & Nepf Reference Tanino and Nepf2008; Xu & Nepf Reference Xu and Nepf2020). In a patchy canopy, the individual patches contribute two length scales; patch width
$W_{p}$
and deflected height
$h$
, from which the patch-wake scale is controlled by
$L_{e}=min (W_{p},h)$
(Marjoribanks et al. Reference Marjoribanks, Lague, Hardy, Boothroyd, Leroux, Mony and Puijalon2019). Specifically, when
$W_{p}\lt h$
, the wake length is dominated by horizontal shear, whereas when
$W_{p}\gt h$
, it is determined by vertical shear. In natural streams, the dimensions of
$W_{p}$
and
$h$
for individual patches have similar scales, as reported in Cornacchia et al. (Reference Cornacchia, Lapetoule, Licci, Basquin and Puijalon2023). This similarity allows either length scale to be used as
$L_{e}$
when analysing the flow around an individual patch. However, in a canopy of multiple patches, if the patches are closely spaced, the turbulent wakes generated by individual patches can merge together downstream from the patches, which can increase the wake length scales (Meire et al. Reference Meire, Kondziolka and Nepf2014). In this case, the wake length is determined by
$h$
. Given this, it is reasonable to assume that, within a patchy canopy,
$L_{e}\sim h$
. However, for patches having a very narrow width (
$W_{p}\ll h$
), this assumption may no longer hold.
Finally, we assume that the flow is fully developed, such that the dissipation rate equals the rate of production (
$\varepsilon =\langle P_{w}\rangle _{z}+\langle P_{s}\rangle _{z}$
). Substituting (2.9), (2.16) and (2.17) into (2.5), the channel-average TKE is
\begin{equation} \left\langle k\right\rangle _{z,pre}={\gamma ^{2}}\left(gS\sqrt{\frac{2gSH}{C_{D}aL_{c}}}+\frac{1.3}{H\sqrt{H\delta _{e}}}K_{c}^{\frac{3}{2}}\big(\tilde{U}_{A}-\tilde{U}_{B}\big)^{3}\left(H-h+\delta _{s}+\delta _{e}\right)\right)^{2/3}h^{2/3}, \end{equation}
in which
$\gamma ^{2}$
is the scale coefficient. Considering that previous studies on emergent vegetation have suggested that
$\gamma ^{2}\sim O(1)$
(
$\gamma ^{2}=1.1\pm 0.2$
, Tanino & Nepf Reference Tanino and Nepf2008;
$\gamma ^{2}=0.96\pm 0.16$
, Xu & Nepf Reference Xu and Nepf2020), this study assumed
$\gamma ^{2}$
to be 1. Note that the variables
$C_{D}$
and
$\delta _{s}$
are expected to be functions of the patch area density, and a comprehensive evaluation of these variables will be described in the result section. Considering that many natural landscapes are characterized by heterogeneous distributions of submerged vegetation patches, rather than homogeneous distributions (see figure 9 in Biggs et al. Reference Biggs, Nikora, Gibbins, Fraser, Green, Papadopoulos and Hicks2018; figure 15 in Biggs, Haddadchi & Hicks Reference Biggs, Haddadchi and Hicks2021; figure 3 in Cameron et al. Reference Cameron, Nikora, Albayrak, Miler, Stewart and Siniscalchi2013; figure 1 in Cornacchia et al. Reference Cornacchia, Van De Koppel, Van Der Wal, Wharton, Puijalon and Bouma2018), the new TKE prediction based on shear and wake turbulence will offer broader applicability compared with previous models that only considered stem-scale wake production in emergent vegetation (Xu & Nepf Reference Xu and Nepf2020) or shear production in a uniform array of submerged stems (Zhang et al. Reference Zhang, Lei, Huai and Nepf2020). However, it is important to note that this model is applicable only to dense canopies, as defined by Belcher, Jerram & Hunt (Reference Belcher, Jerram and Hunt2003), for which the canopy drag is larger than the bed drag, such that a two-layer momentum model is appropriate. Specifically, Belcher et al. (Reference Belcher, Jerram and Hunt2003) defined the dense canopy behaviour for
$C_{D}aL_{c}\geq 0.1$
. For sparse canopies, the shear at the top of the canopy is weak compared with bed shear, so that the velocity profile more closely resembles a turbulent boundary layer, and the two-layer assumption breaks down, leading to an inaccurate prediction of canopy shear production, as defined in this paper. In the present study, all configurations, except for X3Y3 with the smallest
$\phi _{p}$
, satisfied the dense canopy condition.
3. Experimental methods
3.1. Experimental set-up
Laboratory experiments were conducted in a 13-m-long and 0.64-m-wide (
$W$
) recirculating flume. The flow rate was controlled by a variable-speed pump and measured with a flow meter installed in the recirculation pipe (figure 2
a). The water depth was set by a weir located at the downstream end of the flume and measured using rulers attached to the side wall of the flume. For all experiments, the flow rate and the water depth were Q = 450
$\pm$
3 L min–1 and H = 0.193
$\pm$
0.002 m, respectively, such that the channel-averaged velocity
$U_{0}(=Q/HW)$
was 0.061
$\pm$
0.001 m s–1. The test section was located 3.5 m from the flume inlet and had a length (L
test
) that ranged from 2.3 to 4.9 m depending on the spacing between patches.
Figure 2. (a) Schematic diagram of the recirculating flume. Not to scale. Patches were installed on baseboards (grey line) that ran along the length of the test section. The ADV (Vectrino, Nortek) was mounted on a traverse installed above the flume, allowing measurements at desired positions in 3-D space. (b) Side and top view of an enlarged patch. (c) Definition of the gap distance in the x and y directions. (d) A top view of the spatial configuration of each experimental case; yellow characters in the bottom right corner denote the name of the case.
Each patch was constructed from a bundle of 24 flexible plastic plants, Rotala indica, fixed homogeneously to a 5-mm thick PVC plate with area 0.12 m by 0.12 m. Each plant had a stem length of 0.18 m and stem diameter of 0.2 cm. Each stem had 28 leaves distribution along its length. Due to the greater rigidity of model plants compared with real plants (flexural rigidity, EI = 4.2 N m–2 for a model plant, 1.5 N m–2 for a real plant; see table 1 in Zhang & Nepf Reference Zhang and Nepf2020), we artificially deflected the vegetation patches in the direction of flow using a steel wire of 1 mm diameter that wrapped around both the patch and the baseboard (figure 2
b). The deflected patch height, patch width and patch length were
$h$
= 10.6
$\pm$
0.7 cm,
$W_{p}$
= 11.2
$\pm$
0.8 cm and
$L_{p}$
= 26.0
$\pm$
0.9 cm, respectively, which fell within ranges of submerged macrophytes observed by Cornacchia et al. (Reference Cornacchia, Lapetoule, Licci, Basquin and Puijalon2023) in natural rivers (see figure A1 of the supplementary material available at https://doi.org/10.1017/jfm.2025.12). The reconfiguration of patches using steel wire limits the swaying motion that is sometimes present in flexible vegetation, triggered by the passage of turbulent structures (Ghisalberti & Nepf Reference Ghisalberti and Nepf2002; Siniscalchi & Nikora Reference Siniscalchi and Nikora2012). Swaying motion occurs at the top of a canopy (called monami, Ackerman & Okubo Reference Ackerman and Okubo1993) and individual patches may also exhibit swaying motion (Siniscalchi & Nikora Reference Siniscalchi and Nikora2013). The swaying motion can reduce the canopy shear turbulence by weakening the shear layer at the top of the canopy (Nepf Reference Nepf2012). In this regard, canopy shear turbulence generated in our model patches are expected to be larger than that in real vegetation patches, which will be described in the result section.
Table 1. Spatial configurations and flow characteristics. Here
$\sigma _{SE}$
represents the standard error.
Eight spatial configurations were considered, with different longitudinal (
$D_{x}$
) and lateral (
$D_{y}$
) gaps between patches (figure 2
c). The patch locations were determined in two steps. First, the patches were distributed with uniform spacing
$D_{x}$
and
$D_{y}$
. Here,
$D_{x}$
was the distance from the trailing edge of the upstream patch to the upstream edge of the downstream patch (red dashed arrows in figure 2
c), and
$D_{y}$
was the lateral distance between patches (navy dashed arrows in figure 2
c). Second, each patch location was adjusted by a random amount. Specifically, starting from the initial patch location (
$x$
,
$y$
), the adjusted location was (
$x+\alpha _{1}D_{x}$
,
$y+\alpha _{2}D_{y}$
), with
$\alpha _{1}$
and
$\alpha _{2}$
generated using a random number generator in MATLAB, with the constraint
$| \alpha _{1}| ,| \alpha _{2}|\lt 0.3$
to prevent overlap between patches. After the random adjustment, the final gaps
$D_{x}$
and
$D_{y}$
were recalculated, which are described in table 1. To capture both skimming flow and wake interference flow, this study selected gaps in the range of
$0.22\lt D_{x}/h\lt4.83$
and
$0.09\lt D_{y}/W_{p}\lt1.18$
. The cases are denoted by characters ‘X’ and ‘Y’ followed by a nominal number, in which a larger number signifies a larger gap size. We classified the configurations into two types based on lateral spacing between patches: 2-D and 3-D distributions. Specifically, 2-D distributions referred to channel-spanning configurations (X1Y0 and X2Y0) with
$D_{y}/W_{p}\lt0.2$
, while 3-D distributions had
$D_{y}/W_{p}\gt0.2$
. Each test case is illustrated in figure 2(d).
The channel-scale stem density
$m_{stem}$
(
$=n_{p}n_{stem}$
) ranged from 107 to 651 stems m–2, with
$n_{p}$
the number of patches per bed area and
$n_{stem}$
(
${=}24$
) the number of stems in a single patch. The porosity of a single patch (
$0.58=1- \phi _{\textit{single}}$
) was estimated from the solid volume fraction of a single patch (
$\phi _{\textit{single}}$
). The bed area occupied by a single patch was
$A_{p}= 0.028\pm0.001\,{\rm m}^2$
, based on 32 aerial images. The bed area fraction occupied by patches
$\phi _{p} (=n_{p}A_{p})$
was in the range of 0.13–0.79.
3.2. Measurements
Three-dimensional, instantaneous, velocity components, u, v and w, were measured using acoustic Doppler velocimetry (Vectrino, Nortek) (figure 3). The measurements were made across a y–z transect (shaded plane in figure 3 a), because flow adjustment around the patches occurred in both horizontal and vertical planes. The location of the measurement cross-section was chosen to be the region where fully developed flow was achieved. The fully developed flow was determined by confirming that both mean velocity and TKE varied by less than 5 % between different transects along the streamwise direction within the canopies. At every measurement point, the velocity was recorded for 2 min at 100 Hz, which was sufficient to resolve turbulent flow information in the measured data (Park & Hwang Reference Park and Hwang2021). Since ADV measurements are not feasible at less than 5 cm from the water surface (Nortek 2018), the mean and turbulent flow characteristics within this area were obtained by applying linear extrapolation to the measured data (grey areas in vertical profiles in figures 4–6). In post-processing, spikes in the time-series data were eliminated by the method introduced by Goring & Nikora (Reference Goring and Nikora2002) before estimating turbulent quantities.
Figure 3. (a) Schematic of velocity measurements. The yellow shaded plane represents the transect position, and the red circles indicate the measurement resolution in the cross-stream and vertical directions. Relative error in channel-mean velocity and TKE depended on (b) lateral and (c) vertical intervals. The relative error,
$\psi$
, was defined as the ratio of each parameter value,
$\langle k_{t}\rangle _{yz}$
and
$\langle \overline{u}\rangle _{yz}$
, estimated at spacing (
$\Delta y$
,
$\Delta z$
) to the parameter value estimated using the smallest spacing.
Prior to starting the main experiment, a convergence test was performed to determine the optimal spacing between measurement points in the y–z plane. Velocity was measured at tight intervals and used to identify the resolution of measurement at which the transect-averaged streamwise velocity
$\langle \overline{u}\rangle _{yz}$
and TKE
$\langle k_{t}\rangle _{yz}$
became insensitive to the measurement resolution in the spanwise
$(\Delta y)$
and vertical
$(\Delta z)$
directions. The TKE was calculated as
$k_{t}=(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}})/2$
. As shown in figures 3(b) and 3(c), the relative error of
$\langle \overline{u}\rangle _{yz}$
and
$\langle k_{t}\rangle _{yz}$
exceeded
$\pm 5$
% when
$\Delta y/W_{p}$
and
$\Delta z/h$
exceeded 0.4 and 0.2, respectively. Hence, the measurement intervals
$\Delta y=3\,\rm{cm}$
(
$\Delta y/W_{p}=0.27$
) and
$\Delta z=1\,\rm{cm}$
(
$\Delta z/h=0.09$
) were chosen for the remaining transects. This resulted in 19 points in the lateral direction and 14 points in the vertical direction, resulting in a total of 266 points in each transect plane. The number of transects in the gap was determined based on
$D_{x}$
. For “X0Y0”, “X0Y1”, “X0Y2” with small
$D_{x}$
(
$D_{x}/h\approx0.2$
),
$\langle \overline{u}\rangle _{yz}$
and
$\langle k_{t}\rangle _{yz}$
differed by 7.9 % and 9.8 %, respectively, depending on the xposition in the gap; thus, we only measured velocity at a single transect for these cases. For cases with large
$D_{x}/h$
(
$D_{x}/h\gt0.3$
), the velocity transect was carried out at 3–5 different streamwise positions in the gap. The distribution of transects in this manner was reasonable because the turbulent flow characteristics vary linearly from the upstream to downstream within the gap for
$D_{x}/h\lt 5$
(Folkard Reference Folkard2005; Chung et al. Reference Chung, Mandel, Zarama and Koseff2021).
Due to the difficulty in measuring the very small water surface slope (
$S$
) within the test section,
$S$
was calculated from the balance of the momentum flux and the gravitational force in the overflow region (
$h\lt z\lt H$
):
$S=-({\langle \overline{u^{\prime}w^{\prime}}\rangle _{z=h}}/({gH-h}))$
, following King et al. (Reference King, Tinoco and Cowen2012). Here,
$\langle \overline{u^{\prime}w^{\prime}}\rangle _{z=h}$
is the RS at the top of the canopy, which can be calculated from a linear least squares fit to the RS profile above the canopy. Calculation of the water surface slope based on the balance of the momentum flux and the gravitational force in the overflow region is only valid for dense canopies with
$C_{D}aL_{c}\gt0.1$
. The canopy drag coefficient (
$C_{D}$
) was also computed from the momentum equation instead of by direct measurement. Substituting the water surface slope into (2.7), the channel-averaged canopy drag coefficient can be calculated as
\begin{equation}C_{D}=-\frac{2H\big\langle\, \overline{u^{\prime}w^{\prime}}\,\big\rangle _{z=h}}{aL_{c}\left(H-h\right)U_{v}^{2}},\end{equation}
which will be described in the following section.
4. Results
4.1. Evolution of flow structure
Different turbulent flow structures developed within the canopy, depending on the spatial configuration of the 3-D submerged patches. First, we compared the mean velocity and the RS distributions for the configurations X0Y0, X0Y1, X0Y2, which have similarly small longitudinal gap sizes (
$D_{x}/h\lt0.3$
, table 1), but different lateral gap sizes (
$0.1\lt D_{y}/W_{p}\lt1.2$
). For the most dense case, X0Y0, the mean velocity in the vegetated layer (
$0\lt z\lt h$
),
$U_{v}$
, was the lowest (
$U_{v}/U_{0}\lt 0.2$
, figure 4
c) and the velocity within the lateral gaps between patches was similar to that behind the patches (see even blue shading in figure 4
a). Spatial variation in velocity was primarily in the vertical direction, with a strong velocity gradient and high RS at the top of the canopy (
$z/h=1$
, figures 4
c and 4
d). The RS decays sharply towards zero as
$z/h$
approaches 0.5, indicating a limited penetration of vertical turbulent momentum flux into the canopy (figure 4
c). This was the densest canopy considered, and exhibited the smallest penetration distance, and specifically, the mixing layer did not penetrate to the bed. Similar results were observed by Chen et al. (Reference Chen, Jiang and Nepf2013), reporting that dense submerged canopies with
$C_{D}ah\gt0.3$
do not allow the mixing layer to penetrate the full canopy height. The ratio of minimum to maximum velocity (
$\langle \overline{u}\rangle _{min}/\langle \overline{u}\rangle _{max}\approx0.05$
) and the patch area density
$\phi _{p}=0.78$
were consistent with previous criteria for skimming flow. Specifically,
$\langle \overline{u}\rangle _{min}/\langle \overline{u}\rangle _{max}\lt0.05$
(Folkard Reference Folkard2011) and
$\phi _{p}\gt0.4$
(Wolfe & Nickling Reference Wolfe and Nickling1993). Thus, the flow formed in X0Y0 was characterized as skimming flow.
Figure 4. Contours of (a) the time-mean streamwise velocity and (b) the RS within and above the canopy for configurations of X0Y0 (
$\phi _{p}=0.78$
), X0Y1 (
$\phi _{p}=0.47$
) and X0Y2 (
$\phi _{p}=0.32$
). The contours were generated by linearly interpolating 266 measurement points across the y–z transect. Green lines outline the locations of vegetation patches. Vertical distribution of the horizontally averaged (c) streamwise velocity and (d) the RS within the gap. Horizontal bars represent the standard error of the mean across the canopy heterogeneity.
When the lateral gap increased to
$D_{y}/W_{p}=0.63$
(X0Y1), the velocity in the gaps was higher than that in the wake of individual patches (see variation in blue shading in X0Y1, figure 4
a), which was associated with an increase in laterally averaged velocity within the vegetated layer, compared with X0Y0 (compare red and green symbols in figure 4
c). The velocity ratio
$\langle \overline{u}\rangle _{min}/\langle \overline{u}\rangle _{max}=0.19$
indicated a transition away from skimming flow (Folkard Reference Folkard2011). The reduced vertical shear (figure 4
c) resulted in a decrease in RS near the top of the canopy (figure 4
d), compared with the skimming flow case (X0Y0). The small longitudinal gap still restricted vortex penetration into the canopy, resulting in RS approaching zero at
$z/h\approx0.5$
, similar to X0Y0.
When the lateral gap increased further to
$D_{y}/W_{p}=1.18$
(X0Y2), even more flow passed through the gaps, resulting in significant lateral variation of the in-canopy velocity (see variation in blue shading in X0Y2, figure 4
a), as well as a further reduction in the difference between velocity within and above the canopy flow (triangles, figure 4
c), and a further reduction in peak RS (triangles, figure 4
d). Whereas the densest canopy (X0Y0) was associated with a laterally uniform RS (red band in figure 4
b), RS varied laterally in case X0Y2, with significant RS (
$-\langle \overline{u^{\prime}w^{\prime}}\rangle /U_{0}^{2}\gt0.03$
) only observed behind individual patches (figure 4
b).
Figures 4(c) and 4(d) illustrated clear differences in mean velocity and RS distribution depending on patch configuration. An increase in lateral gap
$D_{y}$
(or decrease in
$\phi _{p}$
) decreased the drag discontinuity at the top of the canopy, leading to a reduction in
$\partial \langle \overline{u}\rangle /\partial z$
. This diminished the shear production, which in turn weakened vertical momentum transport. Consequently, the maximum RS value at
$z=h$
decreased with increasing lateral gap length (figure 4
d). Contrary to the region of
$z/h\gt 1$
, the difference in RS profiles becomes negligible as
$z/h$
approaches 0.5, indicating that varying only the lateral gap has little influence on the penetration scale of overflow into the canopy.
Figure 5. Contours of (a) the time-mean streamwise velocity and (b) the RS within and above the canopy for configurations of X1Y1 (
$\phi _{p}=0.45$
,
$D_{x}/h=1.4$
), X2Y2 (
$\phi _{p}=0.28$
,
$D_{x}/h=2.3$
) and X3Y3 (
$\phi _{p}=0.13$
,
$D_{x}/h=4.8$
). The contours were generated by linearly interpolating 266 measurement points across the y–z transect. Green lines outline the locations of vegetation patches. Vertical distribution of the horizontally averaged (c) streamwise velocity and (d) the RS within the gap. The skimming flow case X0Y0 is included for reference. Horizontal bars represent the standard error of the mean across the canopy heterogeneity.
For canopies with 3-D distributed patches, mean and turbulent flow varied systematically with changing
$\phi _{p}$
, which is illustrated by simultaneously varying both longitudinal and lateral gap lengths (figure 5). Contours in figures 5(a) and 5(b) were obtained by averaging across multiple transects within the longitudinal gap. As the gap sizes increased (decreasing
$\phi _{p}$
), more flow passed between patches (figure 5
a), increasing the canopy-average velocity and weakening the vertical shear layer (figure 5
c), which in turn alters the RS (figure 5
d). In particular, as
$\phi _{p}$
decreased, the peak RS decreased in magnitude and shifted downward (figure 5
d). For configuration X3Y3, which has the smallest patch area density,
$\phi _{p}=0.13$
, the canopy drag was too small to generate a shear layer with an inflection point, and the vertical profile of velocity instead resembled a turbulent boundary layer (purple triangles in figure 5
c). This result is more clearly presented in figure A2 of the supplementary material, where the position with the second-order derivative
$\partial ^{2}u/\partial z^{2}=0$
becomes unclear with decreasing
$\phi _{p}$
and eventually disappears for X3Y3. This transition in velocity profile shape for
$\phi _{p}=0.13$
suggested that the boundary, rather than the canopy, was the dominate source of flow resistance in this case, consistent with the description in Belcher et al. (Reference Belcher, Jerram and Hunt2003), who defined a transition from boundary layer flow to canopy flow when the drag contributed by the canopy exceeds the drag contributed by the bed.
Figure 6. Contours of (a) the time-mean streamwise velocity and (b) the RS for channel-spanning 2-D configurations X1Y0 (
$\phi _{p}=0.47$
,
$D_{x}/h=2.3$
) and X2Y0 (
$\phi _{p}=0.32$
,
$D_{x}/h=4.8$
). The contours were generated by linearly interpolating 266 measurement points across the y–z transect. Green lines outline the locations of vegetation patches. (c) Variation in the transect-averaged streamwise velocity in the vegetated layer,
$U_{v}$
, within the gap between channel-spanning patches (see schematic in inset figure). Vertical distribution of the horizontally averaged (d) streamwise velocity and (e) the RS. The 3-D configurations X0Y1 (
$\phi _{p}=0.47$
, light green circle) and X0Y2 (
$\phi _{p}=0.32$
, light purple triangle), which have the same
$\phi _{p}$
as X1Y0 and X2Y0, respectively, were included for comparison. Horizontal bars represent the standard error of the mean across the canopy heterogeneity.
The 2-D channel-spanning configurations X1Y0 and X2Y0 with different longitudinal gaps (
$D_{x}/h=2.3$
and 4.8, respectively) are compared in figure 6. These 2-D configurations produced velocity and RS distributions that were laterally uniform (figures 6
a and 6
b). For non-porous, channel-spanning obstacles, a recirculation zone with negative velocity develops in the wake (Hamed & Chamorro Reference Hamed and Chamorro2018; Le Ribault et al. Reference Le Ribault, Vinkovic and Simoëns2021; You & Tinoco Reference You and Tinoco2023). However, the porous vegetation used in this study allowed sufficient flow through the patches to prevent the formation of a recirculation, i.e. no negative velocity regions were observed (figure 6
a). Within the longitudinal gap, the streamwise velocity in the vegetated layer (
$U_{v}$
) increased monotonically with distance from the upstream patch, as illustrated in figure 6(c).
Comparing figures 4 and 5 (3-D-distributed patches) with figure 6 (2-D patches), it is evident that channel-spanning canopies generate a stronger vertical momentum flux than the unconfined (3-D) canopies. Specifically, while red-coloured regions with
$-\langle \overline{u^{\prime}w^{\prime}}\rangle /U_{0}^{2}\gt$
0.03 were locally measured behind individual patches for 3-D canopies (figures 4
b and 5
b), for channel-spanning patches, this magnitude of RS was present across the entire channel width and penetrated close to the bed (figure 6
b). This resulted in differences in the vertical distributions of horizontally averaged RS. For channel-spanning patches X1Y0 (dark green circle in figure 6
e) and X2Y0 (dark purple triangle in figure 6
e), the maximum RS was close to
$z/h=1$
and decreased linearly to zero at the bed. Compared with a configuration with distributed patches but comparable
$\phi _{p}$
(light green circle and light purple triangle in figure 6
e), the channel-spanning patches exhibited much higher RS within the canopy layer, demonstrating that laterally confined canopies generate stronger vertical momentum transport into the canopy layer.
4.2. Drag and vortex penetration scale in canopy flow
For a canopy of 3-D distributed patches, the canopy drag coefficient increased linearly with increasing
$\phi _{p}$
, ranging from 0.5 to 2.3 (navy, yellow, green dots in figure 7
a), which fell within the typical range of drag coefficient observed in canopies (e.g. Dunn, López & García Reference Dunn, López and García1996; Tang et al. Reference Tang, Tian, Yan and Yuan2014; Park & Hwang Reference Park and Hwang2019; Cui, Felder & Kramer Reference Cui, Felder and Kramer2023). Channel-spanning patches (two dimensional) had higher drag coefficients, compared with 3-D patches with the same
$\phi _{p}$
, which was associated with a stronger sheltering effect between distributed patches. Sheltering occurs when a downstream element is situated within the wake created by an upstream element, such that the downstream element is exposed to a reduced velocity and, thus, reduced drag (Etminan, Lowe & Ghisalberti Reference Etminan, Lowe and Ghisalberti2017). The sheltering effect varies with element arrangement and is the most pronounced for in-line configurations (Sumner Reference Sumner2010). Consequently, drag reduction is primarily achieved by decreasing the longitudinal gap size rather than the lateral gap size. Specific to this study, for the same
$\phi _{p}$
, channel-spanning patches had the largest longitudinal gap length, and thus, experienced the least sheltering effect, resulting in the highest canopy drag coefficient (red dots in figure 7
a).
Figure 7. The variation in (a) the drag coefficient
$C_{D}$
and (b) the position of maximum RS
$\tilde{h}$
versus patch area density
$\phi _{p}$
. (c) Comparison between predicted penetration length scale (
$\delta _{e,pre}$
, (4.2)) and observed
$\delta _{e,obs}$
based on the transect-averaged vertical RS profiles in figures 4–
6. Vertical bars in (a) indicate the uncertainty in
$C_{D}=-{2H\langle \overline{u^{\prime}w^{\prime}}\rangle _{z=h}}/{aL_{c}(H-h)U_{v}^{2}}$
. Horizontal bars in (c) represent the uncertainty in (4.2) and vertical bars indicate the uncertainty in measurement.
The variation in canopy drag coefficient due to spatial arrangement diminished as
$\phi _{p}$
increased. Specifically, the maximum difference in
$C_{D}$
was 1.2 at
$\phi _{p}=0.3$
between X0Y2 and X2Y0, but only 0.6 at
$\phi _{p}=0.5$
between X0Y1 and X1Y0 (figure 7
a). At lower
$\phi _{p}$
, the sheltering effect between patches is more sensitive to patch arrangement, because a greater range of spacing is possible, leading to greater variation in individual patch drag. However, as
$\phi _{p}$
increases, more of the patches experience a similar degree of sheltering, because the range of spacing is more limited. This result was similar to Gijón Mancheño et al. (Reference Gijón Mancheño, Jansen, Winterwerp and Uijttewaal2021), who demonstrated that cylinder arrays with higher porosity exhibited larger variation in bulk drag coefficient due to arrangement compared with those with lower porosity. This supports that the role of sheltering in producing variations in drag is similar across different object morphologies, producing a greater sensitivity to object arrangement at higher canopy porosity. Eventually,
$\phi _{p}$
is high enough to produce skimming flow, for which every patch was sheltered within the vegetation layer, and configuration had no impact on
$C_{D}$
. In summary, the significance of patch arrangement to flow structure and canopy drag weakens with increasing
$\phi _{p}$
, eventually becoming negligible for skimming flow.
A linear relation between canopy drag coefficient and
$\phi _{p}$
can be approximately represented by
$C_{D}=(3.0\pm 0.7)\phi _{p}$
(with a 95 % confidence interval (2.3, 3.7)). This relation is valid for patch Reynolds number,
$Re_{p}$
(
$=U_{0}h/\upsilon$
, in which
$\upsilon$
is kinematic viscosity) equal to 6000, but the slope of the linear function in figure 7(a) should depend on Reynolds number, and thus on the submergence ratio (
$H/h$
). Specifically, as the submergence ratio increases, the velocity (and patch Reynolds number) within the canopy decreases, leading to earlier flow separation (more upstream) around the individual patches and, thus, wider wakes. As a result, as
$U_{v}$
and thus Reynolds number decreases, the vegetation patches experience higher drag, which is reflected in a higher drag coefficient. This means that as
$U_{v}$
decreases with increasing
$H/h$
, the drag coefficient
$C_{D}$
increases. Consequently, the slope of the linear function between
$C_{D}$
and
$\phi _{p}$
was expected to be steeper with increasing
$H/h$
or decreasing the Reynolds number.
The maximum RS occurred at a distance
$\tilde{h}$
from the bed (figure 1). From figures 4 to 6, the location of the maximum RS decreased from
$\tilde{h}/h\approx 1\ \textrm{to}\ 0.7$
with decreasing patch density. If we additionally include the limits of the bare bed (
$\phi _{p}=0$
), for which the maximum RS is at the bed (
$\tilde{h}/h=0$
) and a continuous canopy (
$\phi _{p}=1$
), for which
$\tilde{h}/h=1$
, we find the relationship
$\tilde{h}/h=\phi _{p}^{0.11\pm 0.02}$
(95 % CI, and R-squared
$=$
0.99, figure 7
a), which will be used to predict shear production in § 4.4. The relationships presented in figures 7(a) and 7(b) do not capture the effects of the spatial arrangement of patches as the coefficients were derived solely as a function
$\phi _{p}$
.
Because the patch distribution influences the degree of sheltering and, thus, the total drag provided by the canopy layer, it also influences the penetration of the shear-layer turbulence into the canopy. Nepf et al. (Reference Nepf, Ghisalberti, White and Murphy2007) derived an inverse relation between the non-dimensional drag scale and the penetration scale as
$\delta _{e}/h\sim [C_{D}ah]^{-1}$
. Due to the gap between patches, in the 3-D submerged canopies, the vortex scale shifted down from the top of the canopy by length scale
$\delta _{s}=h-\tilde{h}$
(figure 1
c), and the deflected part of the canopy contributed additional drag. Considering these modifications, we propose that the vortex penetration scale in a canopy of 3-D submerged patches can be described by
\begin{equation}\frac{\delta _{e}}{h}=\frac{0.3\pm 0.1}{C_{D}aL_{c}}-\frac{\delta _{s}}{h}=\frac{0.3\pm 0.1}{C_{D}ah\left(1+\cos ^{2}\theta \right)}-\frac{\delta _{s}}{h}.\end{equation}
Here, the constant scale
$0.3\pm0.1$
was based on Ghisalberti & Nepf (Reference Ghisalberti and Nepf2009). Considering that the penetration scale is additionally constrained by the canopy height
$\tilde{h}$
, in this study, with
$\theta =22.5^{\circ}$
,
\begin{equation} \delta _{e}=\textit{min} \left(\frac{0.17\pm 0.05}{C_{D}a}-\delta _{s},\tilde{h}\right), \end{equation}
in which
$\delta _{s}$
was computed as
$\delta _{s}=h-\tilde{h}=h(1-\phi _{p}^{0.11\pm 0.02})$
, based on the relationship in figure 7(b). The comparison between (4.2) and the observed values is depicted in figure 7(c). The predicted
$\delta _{e}$
have a good agreement with observed values with the mean difference ratio of 21 %. Since the prediction of
$\delta _{e}$
was developed based on a continuum model, it performed better for higher
$\phi _{p}$
compared with lower
$\phi _{p}$
.
4.3. Turbulence production
Within a submerged patchy canopy, turbulence may be generated at various scales. To characterize the turbulent length scale within the canopy, we computed the integral length scales for each measurement position by multiplying the temporal-averaged velocity by the integral time scale and determined
$L_{t}$
by horizontally averaging them as
$L_{t}=\langle \overline{u}T_{i}\rangle$
, in which
$T_{i}$
was the integral time scale obtained from the autocorrelation function following O’Neill et al. (Reference O’Neill, Nicolaides, Honnery and Soria2004). For all experimental cases, the integral length scales within the canopy (
$z/h\lt 1$
) were smaller than those above the canopy (
$z/h\gt 1$
) because the gap between patches confined the development of large-scale turbulence (see vertical profiles of
$L_{t}$
in figures 8
a, 8
b and 8
c). As the patch distribution became more sparse (decreasing patch area density), the increased gap size between patches allowed larger turbulent motions within the canopy. For example, the configuration X0Y0 with the highest
$\phi _{p}$
exhibited the smallest
$L_{t}$
, while X3Y3 with the lowest
$\phi _{p}$
exhibited the largest
$L_{t}$
(red and purple markers, respectively, in figure 8
b). The channel-averaged integral length scale ranged from 2 to 5 cm across different configurations. This range of scales could be contributed from either the leaves (approximately 2 cm) or patch (approximately 10 cm), but contributions from the stem (approximately 2 mm) were unlikely. However, the integral length scale obtained from a single-point velocity measurement was insufficient to determine which element scale (leaves or patch) governed turbulence production within the patchy canopy. Thus, this study did not specify a dominant turbulence scale.
Figure 8. Vertical distribution of the integral length scale: (a,b) 3-D distributed patches and (c) 2-D channel-spanning patches. Variation of channel-averaged shear (blue), wake (red) and total (orange) turbulence production with patch area density. (d,e) Three-dimensional distributed patches and (f) 2-D channel-spanning patches. Horizontal bars in (a–c) represent the standard error of the mean across the canopy heterogeneity. Vertical bars in (d–f) indicate the uncertainty of production terms.
Figure 9. (a) Normalized channel-averaged TKE (
$\langle k_{t}\rangle _{z}/U_{0}^{2}$
) and (b) the relative contribution of shear (
$\Omega _{s}$
, triangle) and wake production (
$\Omega _{W}$
, square) to total turbulence production with respect to
$\phi _{p}$
.
Canopy turbulence was generated by two sources: shear and wake production (figure 1), whose contributions vary depending on the spatial heterogeneity. Figure 8(d–f) compares the channel-averaged shear and wake production, calculated from (2.3) and (2.6), respectively. For configurations with 3-D patches (figures 8
d and 8
e), with increasing
$\phi _{p}$
both the shear and RS at the top of the canopy increased (figures 4 and 5), which increased the channel-averaged shear production (blue markers, figures 8
d and 8
e). With increasing
$\phi _{p}$
, canopy drag (
$C_{D}a$
) increased but mean velocity within the vegetated layer (
$U_{v}$
) decreased. These opposing effects limited the increase in wake production for
$\phi _{p}\gt0.5$
, as depicted by red markers.
The configurations with channel-spanning patches, X1Y0 and X2Y0, exhibited a different trend from the other cases (figure 8
f). Specifically, compared with canopies of 3-D distributed patches with the same
$\phi _{p}$
, canopies consisting of channel-spanning patches had higher turbulence production. In addition, the channel-spanning patches with
$\phi _{p}=0.47$
had higher production than the highest density of distributed canopy (X0Y0,
$\phi _{p}=0.78$
). This was because the significant decrease in
$U_{v}$
associated with skimming flow at
$\phi _{p}=0.78$
decreased wake production (red markers in figure 8
f). Therefore, considering X0Y0 to be a limiting case of channel-spanning patches, for channel-spanning patches, the peak in turbulence production was not associated with the peak vegetation density (figure 8
f).
For cases considered in this study, the channel-averaged TKE increased with
$\phi _{p}$
(figure 9
a). The channel-averaged TKE at
$\phi _{p}=0$
(
$\langle k_{t}\rangle _{z}/U_{0}^{2}=0.01$
, grey circle in figure 9
a) was computed from measurements over a bare bed and was consistent with previous bare bed studies (Nakagawa, Nezu & Ueda Reference Nakagawa, Nezu and Ueda1975; Nezu Reference Nezu1977). The limit of a channel fully covered by patches (
$\phi _{p}=1$
) was also evaluated as follows. Assuming a high solid volume fraction within the individual patches, the velocity within the canopy would be close to zero, resulting in negligible wake production, and thus, turbulence would only be produced by shear at the top of the canopy. With zero velocity in the canopy, we can draw an analogue between open channel flow and the flow above the canopy, such that we expect the turbulence to scale on the velocity above the canopy,
$U_{A}$
(figure 1), then the channel-averaged TKE at
$\phi _{p}=1$
can be estimated to be
$\langle k_{t}\rangle _{A}/{U_{A}}^{2}=0.01$
. Given the submergence ratio
$H/h=1.82$
,
$U_{A}$
was estimated from continuity as
$U_{A}=({1.82}/{0.82})U_{0}\approx 2.2 U_0$
, thus, the channel-averaged TKE at
$\phi _{p}=1$
was approximated as
$\langle k_{t}\rangle _{z}/U_{0}^{2}\approx 0.02$
(purple marker in figure 9
a). This result suggests that the variation of
$\langle k_{t}\rangle _{z}/U_{0}^{2}$
with
$\phi _{p}$
follows a downward concave trend (grey-dotted line in figure 9
a). This physically implies that upon approaching
$\phi _{p}=1$
, the velocity within the canopy becomes low and uniformly distributed, shifting the dominant source of wake production from patch-scale turbulence to individual plant-scale turbulence, which eventually results in lower TKE. For the same
$\phi _{p}$
, decreasing the lateral gap length between patches, tending toward channel-spanning patches, led to stronger turbulence in the canopy (compare navy to orange to red dots in figure 9
a), which was consistent with higher production measured for channel-spanning patches (figure 8).
The relative contribution of wake
$(\Omega _{W}=\langle P_{w}\rangle _{z}/(\langle P_{w}\rangle _{z}+\langle P_{s}\rangle _{z}))$
and shear production
$(\Omega _{S}=\langle P_{s}\rangle _{z}/(\langle P_{w}\rangle _{z}+\langle P_{s}\rangle _{z}))$
to turbulence generation varied with
$\phi _{p}$
. Specifically, at lower
$\phi _{p}$
, wake production (square markers in figure 9
b) contributed more to turbulence production than shear production (triangle markers in figure 9
b), but this trend reversed at higher
$\phi _{p}$
. In addition, the peak TKE was observed at
$\phi _{p}\approx 0.5$
(figure 9
a), for which shear and wake production made nearly equal contributions (
$\Omega _{s}=\Omega_{W}=0.5$
, figure 9
b).
Since this study constrained the movement of individual plants using steel wire, the shear production in our model canopy can differ from that in real submerged canopies with the swaying motion. Specifically, Ghisalberti & Nepf (Reference Ghisalberti and Nepf2006) compared the RS profile for the same flexible submerged canopy under conditions with and without swaying motion and observed a reduction of peak RS by 38 %. Considering the relationship between shear production and the peak RS,
$P_{s,z=\tilde{h}}\sim (-\overline{u^{\prime}w^{\prime}}_{z=\tilde{h}})^{3/2}$
(see (2.10), (2.13) and (2.14)), a 38 % decrease in peak RS could result in a 51 % reduction in the shear production. In this study, the relative contribution of shear production to turbulence production ranged from 30 (
$\phi _{p}=0.13$
) to 100 % (
$\phi _{p}=1$
) (figure 9
b), indicating that the total turbulence production in our model patches might be 15–51 % larger than that in real submerged patches.
4.4. Prediction of channel-averaged TKE
The prediction of channel-averaged TKE must begin with predictions of wake and shear production. Thus, we first verified the prediction of wake (2.9) and shear (2.16) production with measurements, as depicted in figure 10. Predicted wake production agreed with measurements within a relative difference of 9.4 %, except for the configuration X2Y0 (figure 10
a), for which the prediction overestimated the measured value by 37 %. This was because a significant underestimation of the drag coefficient (see figure 7
a) led to an overestimation of
$U_{v}$
.
Figure 10. Comparison of predicted and measured channel-averaged (a) wake production and (b) layer velocities
$\tilde{U}_{A}$
and
$\tilde{U}_{B}$
. Inset of (b) represents
$U_{\tilde{h}}$
for all experimental cases. (c) Predicted
$U_{\tilde{h}}$
for X0Y0 depending on submergence ratio. The data in the inset of (c) were obtained from table 1 in Chen et al. (2013). (d) Variation in TKE depending on
$K_{c}$
. The x and y axes in (a), (b) and (d) indicate the predicted and measured values, respectively. The black dash–dotted lines in (a), (b) and (d) represent
$y=x$
lines.
To predict shear production, the velocity above (
$\tilde{U}_{A}$
, (2.15)) and below (
$\tilde{U}_{B}$
, (2.12))
$z=\tilde{h}$
must be predicted accurately. This study assumed
$U_{\tilde{h}}\approx U_{0}$
to predict
$\tilde{U}_{B}$
from (2.12), because measured
$U_{\tilde{h}}$
values were within 10 % of
$U_{0}$
for all experimental cases (see inset of figure 10
b). The velocity predictions exhibited a good agreement with observed values, with mean deviations of 3.7 % and 9.7 %, respectively (figure 10
b).
Note that the simplifying assumption
$U_{\tilde{h}}\approx U_{0}$
was specific to the canopy submergence
$H/h\approx 2$
, and will likely vary with variation in the submergence ratio. Specifically, data from Chen et al. (Reference Chen, Jiang and Nepf2013) showed that the velocity within the canopy (
$U_{v}$
, see figure 1) exhibited a power-law dependence on the submergence ratio
${U_{v}}/{U_{0}}=({H}/{h})^{q}$
, in which the exponent
$q$
depended on the canopy configuration (figure 10
c–1). For a dense canopy with
$\tilde{h}/h\approx 1$
(figure 7
b),
$\tilde{U}_{B}\approx U_{v}$
such that
$\tilde{U}_{A}$
can be estimated as
$\tilde{U}_{A}\approx U_{A}=(({U_{0}H}/({H-h})(1-({H}/{h})^{q-1})$
from (2.15). Combining (2.10), (2.11) and (2.12),
$U_{\tilde{h}}\approx \tilde{U}_{B}+2.6\sqrt{K_{c}({\delta _{e}}/{H})^{1/3}}(\tilde{U}_{A}-\tilde{U}_{B})$
, and substituting
$\tilde{U}_{A}$
and
$\tilde{U}_{B}$
into this equation,
$U_{\tilde{h}}$
for a dense canopy can be described as a function of submergence:
\begin{equation} \frac{U_{\tilde{h}}}{U_{0}}\approx \left(\frac{H}{h}\right)^{q}+2.6K_{c}^{1/2}\left(\frac{\delta _{e}}{H}\right)^{1/6}\left[\frac{H}{H-h}\left(1-\left(\frac{H}{h}\right)^{q-1}\right)-\left(\frac{H}{h}\right)^{q}\right]. \end{equation}
For the configuration X0Y0 satisfying the assumption of
$\tilde{h}/h\approx 1$
, the value of
$q = -1.73$
(navy dash–dotted line in figure 10
c), and the variation of normalized
$U_{\tilde{h}}$
with submergence ratio is depicted as an orange dash–dotted line in figure 10(c). As submergence increased,
$U_{\tilde{h}}$
decreased while passing
$U_{\tilde{h}}/U_{0}=1$
at
$H/h=2$
, and eventually for terrestrial canopies with
$H/h=\infty$
,
$U_{\tilde{h}}$
will tend to zero in (4.3).
The predicted shear production was also sensitive to the empirical constant
$K_{c}$
(see (2.16)), which can depend on the vegetation characteristic. Specifically, previous studies by Chen et al. (Reference Chen, Jiang and Nepf2013) and Zhang et al. (Reference Zhang, Lei, Huai and Nepf2020) suggested
$K_{c}=0.07\pm 0.02$
for rigid cylinders and
$K_{c}=0.05\pm 0.01$
for flexible plants. In this study, the predicted shear production agreed best with measured values using the rigid canopy value,
$K_{c}=0.07\pm 0.02$
(dark green markers in figure 10
d), which was consistent with the fact that the bundle of stems had negligible movement and, thus, behaved as rigid stems. This result highlights a limitation of our model patches in representing real vegetation patches that exhibit swaying motion. If the motion of individual plants were unrestricted, allowing swaying motion to occur in our model patches, the observed shear production would decrease, potentially leading to better agreement between the predictions and observations for smaller
$K_{c}.$
In addition, the predicted shear production underestimated the measured values, which can be attributed to the exclusion of bed-generated turbulence in our prediction model.
Based on the predicted wake and shear production, the channel-averaged TKE was estimated using (2.18) and compared with the measured values in figure 11(a). The predicted TKE agreed best with the measured values when
$\gamma ^{2}=0.91\pm0.08$
, which was within 10 % of
$\gamma ^{2}=1$
. This coefficient was the same within uncertainty as the values of
$\gamma ^{2}$
in previous studies using rigid emergent vegetation (
$\gamma ^{2}=1.1\pm 0.2$
Tanino & Nepf Reference Tanino and Nepf2008) and emergent flexible vegetation (
$\gamma ^{2}=0.96\pm 0.16$
Xu & Nepf Reference Xu and Nepf2020). This indicates that our model accurately captured the wake and shear turbulence generated by a canopy with a complex 3-D morphology. In addition, the consistent value of
$\gamma ^{2}\approx 1$
, regardless of vegetation type, demonstrated the applicability of our model to diverse types of plants without additional fitting between predictions and measurements.
Figure 11. (a) Normalized predicted (x axis) and observed (y axis) channel-averaged TKE (
$\langle k_{t}\rangle _{z}/U_{0}^{2}$
). The black solid line represents the linear best fit
$\gamma ^{2}=0.91\pm 0.08$
. (b) The variation in channel-average TKE with submergence ratio, based on case X0Y0. Here TKEWake (red line) and TKEShear (blue line) are computed from only the wake production term (thus
$\langle P_{s}\rangle _{z}=0$
) and the shear production term (thus
$\langle P_{w}\rangle _{z}=0$
), respectively. Vertical bars indicate the standard error of measured TKE. Horizontal bars represent the uncertainty in (2.18).
In natural river systems, variability in flow conditions or vegetation growth can change the canopy submergence ratio, which will then alter turbulence level. Thus, we evaluated the applicability of our model to different submergence ratios. For simplicity, we assumed a dense canopy with peak RS at the top of the canopy (
$\delta _{s}=0$
), for which
$\sqrt{{\delta _{e}}/{H}}K_{c}^{({3}/{2})}(\tilde{U}_{A}-\tilde{U}_{B})^{3}\approx u_{*,z=h}^{3}$
(see (2.10) and (2.11)). Substituting
$S={(u_{*,z=h}^{2})}/$
$({g(H-h})$
and
$L_{c}=h(1+\cos ^{2}\theta )$
into (2.9) and (2.18), shear and wake production can be simplified as
\begin{equation} \left\langle P_{s}\right\rangle _{z}=\frac{1.3u_{*,z=h}^{3}}{H\delta _{e}}\left(H-h+\delta _{e}\right), \end{equation}
\begin{equation} \left\langle P_{w}\right\rangle _{z}=\frac{u_{*,z=h}^{3}}{\left(H-h\right)^{1.5}}\frac{H}{h\beta }\sqrt{\frac{2}{\alpha }}.\qquad \end{equation}
Here,
$\alpha =n_{\textit{total}}C_{D}W_{p}/L_{test}W=\phi _{p}C_{D}/L_{p}$
(
$\rm m^{-1}$
) and
$\beta =(1+\cos ^{2}\theta )=1.85$
in this study. Using the velocity at the top of the canopy,
$U_{\tilde{h}}$
, described as a function of submergence ((4.2)) to predict the shear velocity at the top of the canopy ((4.4) and (4.5)), (2.18) can describe the channel-averaged TKE as a function of submergence ratio:
\begin{align} \frac{\left\langle k\right\rangle _{z,pre}}{U_{0}^{2}} & = {\gamma ^{2}}{K_{c}}{\left(\frac{\delta _{e}}{H}\right)^{\frac{1}{3}}}{\left[\frac{H}{H-h}\left(1-\left(\frac{H}{h}\right)^{q-1}\right)-\left(\frac{H}{h}\right)^{q}\right]^{2}}\nonumber\\& \quad \times \left(\frac{1.3}{H\delta _{e}}\left(H-h+\delta _{e}\right) +\frac{1}{\left(H-h\right)^{1.5}}\frac{H}{h\beta }\sqrt{\frac{2}{\alpha }}\right)^{\frac{2}{3}}h^{\frac{2}{3}} . \end{align}
As an example, configuration X0Y0 satisfies the assumption
$\delta _{s}=0$
and, for this case,
$q=-1.73$
(figure 10
c), and
$\alpha$
and
$\delta _{e}$
are 6.9 and 0.33
$h$
, respectively. Substituting these values into (4.6), the change in TKE with submergence ratio
$H/h$
is shown in figure 11(b).
For emergent vegetation (
$H/h=1$
), the contribution of shear production (blue dotted line in figure 11
b) was negligible compared with that of wake production (red solid line in figure 11
b), such that turbulence production within the canopy was governed by wake production. With increasing submergence ratio, the reduction in mean velocity within the canopy led to a sharp decrease in wake production, resulting in comparable wake and shear productions at
$H/h\approx2$
. As submergence increased further (
$H/h\gt3$
), TKEShear+Wake exhibited a similar trend to TKEShear, signifying negligible wake production. Eventually, for terrestrial canopies with
$H/h=\infty$
, wake production in (4.6) approached zero, leading shear production as the sole source of turbulence production.
To evaluate the applicability of our model to different canopy conditions, we used data from Zhao & Nepf (Reference Zhao and Nepf2024), who conducted experiments with a submerged canopy of rigid model vegetation (vertical,
$\theta =90^{\circ}$
, circular cylinders) under different densities (frontal area
$a$
= 1.1 and 4.9
${\rm m}^{-1}$
), submergence ratios (1.4 < H/h< 3.5) and flow conditions (0.33 <
$U_{0}$
<0.4 m s–1). Similar to the TKE model trend in figure 11(b), Zhao & Nepf (Reference Zhao and Nepf2024) also observed that channel-averaged TKE decreased with increasing submergence ratio. To apply our model to this different vegetation morphology, we adjusted the three parameters (
$\alpha$
,
$\beta$
and
$q$
) in (4.6). Specifically,
$\alpha$
and
$\beta$
were calculated as
$\alpha =n_{\textit{total}}C_{D}W_{p}/L_{test}W=C_{D}aH/h$
and
$\beta =(1+\cos ^{2}\theta )=1$
, respectively. Note that, as vegetation density increased, the velocity within the canopy decreased more rapidly with increasing submergence ratio (figure 10
c–1), leading to a decrease in
$q$
. Specifically,
$q=-0.08a-0.27$
, based on data in figure 10(c–1). This fit was used to estimate
$q$
for a = 1.1 and 4.9
${\rm m}^{-1}$
, respectively. The predicted and observed TKE had good agreement (black line in figure 11
a), demonstrating the extension of our model to different vegetation types and submergence ratios. However, the linear relationship between
$q$
and
$a$
was validated only for homogeneously distributed rigid canopies, and
$q$
may need to be recalibrated for other vegetation types. In addition, while we validated our model with different submergence ratios using the data from Zhao & Nepf (Reference Zhao and Nepf2024), this validation has limitations in fully evaluating the model applicability to the diverse real-world conditions with 3-D submerged patchy canopies. Validation under diverse conditions of submerged patchy canopies would be beneficial to ensure the broader applicability of our model.
Many physical processes in natural catchments are influenced by vegetation-generated turbulence, so that it is crucial to accurately model and measure turbulence levels in the field. The development of aerial photography tools such as drones and satellite imaging have facilitated the acquisition of the large-scale spatial distribution of vegetation in natural streams. These observations could be combined with the prediction described in this study to rapidly assess turbulence levels without requiring velocity measurements and over the scale of the entire river network. Furthermore, the finding in figure 11(b) can support simplifications in TKE modelling, depending on the submergence ratio. Specifically, when canopies in the test area are close to emergent, the model can be simplified by neglecting shear production, whereas for deeply submerged canopies, neglecting wake production is reasonable.
Using patch area density
$\phi _{p}$
as a representative parameter in our model could be effective to predict TKE in vegetated channels based on aerial images, but this approach contains limitations that depend on plant morphology. Specifically, for plants with a vertically non-uniform distribution of frontal area, such as P. natans and Odorata, patch area density might be overestimated due to their wide leaves. This plant morphology will create weaker vegetation drag and a less pronounced shear layer, leading to weaker turbulence and potentially overestimating TKE. In this regard, the applicability of our model is limited to plants with a vertically uniform distribution of frontal area, such as Callitriche platycarpa (see figure 1 in Cornacchia et al. Reference Cornacchia, Van De Koppel, Van Der Wal, Wharton, Puijalon and Bouma2018, 2019), G. densa (Cornacchia et al. Reference Cornacchia, Riviere, Soundar Jerome, Doppler, Vallier and Puijalon2022). In addition, it should be noted that our model is only applicable for dense canopies, satisfying a two-layer model, where there is a significant difference in velocity between within and above the canopy. The accuracy of our model may decrease for sparse vegetation because the flow evolution is predominantly governed by the bottom boundary layer rather than the canopy shear layer.
5. Conclusion
Turbulence production within a patchy submerged canopy was investigated under a range of spatial configurations with different longitudinal and lateral gap sizes. With increasing patch area density (decreasing distance between patches), the difference in mean velocity within and above the canopy increased, enhancing the shear layer at the top of the canopy. The decrease in velocity within the canopy and enhancement of shear layer shifts the main source of turbulence production from wake production to shear production. Specifically, at low patch density with relatively high mean velocity within the canopy and a weak shear layer, wake production was the dominant source in turbulence production compared with shear production, whereas for a high patch density, shear production became dominant. For the same patch area density, 2-D channel-spanning patches exhibited deeper turbulence penetration into the canopy and generated larger turbulence compared with 3-D-distributed patches. A new prediction model for the channel-averaged TKE provided a good prediction for the submerged patchy canopies with heterogeneous distribution (figure 11 a). Finally, the turbulent flow characteristics observed in this study were also similar to those observed around coral/mussel mounds (e.g. Constantinescu, Miyawaki & Liao Reference Constantinescu, Miyawaki and Liao2013; Chen et al. Reference Chen, He, Fang, Liu and Dey2024). In this regard, we anticipate that our findings can be applied to enhance the understanding of flow dynamics in coral/mussel mound communities.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jfm.2025.12.
Funding
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2022R1A6A3A03068939).
Declaration of interests
The authors report no conflict of interest.
