1. Introduction
Canopy flows develop over impermeable surfaces where slender, attached and often flexible elements form a complex interfacial layer. Such flows are ubiquitous in both natural and built environments, including terrestrial and aquatic vegetation, as well as urban landscapes. In forests, canopy–flow interactions regulate the exchange of carbon dioxide and oxygen (Bui et al. Reference Bui, Wallace, Kavassalis, Alwe, Flynn, Erickson, Alvarez, Millet, Steiner and Griffin2021); in aquatic systems, vegetation supports biodiversity and mitigates coastal hazards by damping wave energy (Lei & Nepf Reference Lei and Nepf2019). In urban contexts, understanding canopy-layer dynamics can inform building designs that improve pollutant dispersion (Michioka et al. Reference Michioka, Takimoto and Sato2014).
Canopies are typically classified as emergent (
$h/H\gt 1$
) or submerged (
$h/H\lt 1$
), depending on the ratio between the element height
$h$
and the fluid depth
$H$
. In emergent canopies, the flow is dominated by filament-wake eddies throughout the entire fluid depth (van Rooijen et al. Reference van Rooijen, Lowe, Ghisalberti, Conde-Frias and Tan2018). In contrast, submerged canopies exhibit multiscale dynamics, characterised by the coexistence of large-scale motions and smaller-scale turbulence (Monti et al. Reference Monti, Nicholas, Omidyeganeh, Pinelli and Rosti2022). For submerged configurations, fluid–structure interactions are strongly influenced by the solidity
$\lambda$
, i.e. the ratio of frontal to bed area occupied by the canopy (Nepf Reference Nepf2012). In dense canopies (high
$\lambda$
), closely spaced elements disrupt the coherence of incoming eddies, limiting their penetration (Sharma & García-Mayoral Reference Sharma and García-Mayoral2020b
), whereas in sparse canopies (low
$\lambda$
) larger spacings allow overlying turbulent motions to extend down to the wall (Sharma & García-Mayoral Reference Sharma and García-Mayoral2020a
). Recent findings suggest that solidity alone is insufficient to characterise canopy-induced dynamics, and that additional geometrical parameters must be considered (Monti et al. Reference Monti, Nicholas, Omidyeganeh, Pinelli and Rosti2022). In this study, we focus on dense submerged canopies composed of slender filaments (
$h \gg d$
, with
$d$
the cross-sectional stem scale), representative of aquatic vegetation.
Turbulent flow over dense, submerged canopies exhibits pronounced vertical stratification (Poggi et al. Reference Poggi, Katul and Albertson2004). Analogously to flows over porous media, the domain can be broadly divided into two regions: an inner region within the canopy, where velocity fluctuations are attenuated, and an outer region above the canopy, where the mean velocity profile approaches a logarithmic-like behaviour (Chen & García-Mayoral Reference Chen and García-Mayoral2023). At the canopy tips, a shear layer develops due to the abrupt change in drag (Raupach et al. Reference Raupach, Finnigan and Brunei1996), which is susceptible to Kelvin–Helmholtz (KH)-like instabilities (Finnigan Reference Finnigan2000). The resulting spanwise-coherent KH rolls evolve under mean shear into more complex, three-dimensional hairpin-like vortices, which eventually break down into smaller turbulent structures (Bailey & Stoll Reference Bailey and Stoll2016). These large-scale motions coexist with smaller-scale turbulence and extend their influence into the canopy interior (Sharma & García-Mayoral Reference Sharma and García-Mayoral2020a ; Monti et al. Reference Monti, Nicholas, Omidyeganeh, Pinelli and Rosti2022; Foggi Rota et al. Reference Foggi Rota, Monti, Olivieri and Rosti2024). Observations show that the penetration and coherence of these motions are affected by element spacing, stem inclination and flexibility.
Over the past few decades, numerous studies have investigated turbulence in canopy flows using single-point turbulent kinetic energy budgets (Wilson & Shaw Reference Wilson and Shaw1977; Shen & Leclerc Reference Shen and Leclerc1997; López & García Reference López and García2001; Poggi et al. Reference Poggi, Katul and Albertson2004; Defina & Bixio Reference Defina and Bixio2005; Cui & Neary Reference Cui and Neary2008). These works show that turbulent kinetic energy is primarily generated near the canopy tip, where mean shear is strongest, and subsequently redistributed both within the canopy and into the overlying flow. However, single-point measurements cannot capture the scale dependence of these processes or identify the flow structures responsible for production, redistribution and inter-scale transfer. More recently, He et al. (Reference He, Liu and Shen2022) applied a spectral turbulent kinetic energy budget to relate energy distribution to streamwise Fourier modes. While this provides a first step towards a scale-aware perspective, it remains limited to one-dimensional modes and turbulent kinetic energy alone, leaving three-dimensional interactions, anisotropy and connections to coherent structures unresolved.
Despite evidence of multiscale dynamics, the drivers of turbulence in dense canopies remain unclear. In particular, the scales dominating production, pressure–strain redistribution and energy exchange across the canopy-tip interface are not yet quantified. A detailed, scale-resolved characterisation connecting energy dynamics to coherent structures is therefore needed. The objective of this work is indeed to identify the scale-dependent mechanisms that sustain turbulence in dense submerged canopies, and to determine how energy is exchanged between the canopy layer and the overlying flow. Building on the direct numerical simulation (DNS) database of Monti et al. (Reference Monti, Olivieri and Rosti2023), we analyse the Reynolds stresses budget in the joint space of wall-normal position and separation scale. We employ the anisotropic generalised Kolmogorov equations (AGKE) (Gatti et al. Reference Gatti, Chiarini, Cimarelli and Quadrio2020), which provide scale-space transport equations for each Reynolds stress component and enable a local characterisation of production, pressure–strain redistribution and inter-scale fluxes. This framework allows us to determine which scales are dynamically active in different regions of the flow, how they mediate canopy–outer coupling and how these mechanisms relate to coherent structures. Extending the spectral analysis of He et al. (Reference He, Liu and Shen2022), we resolve the full tensorial Reynolds stresses budget and consider both streamwise and spanwise separations, providing a three-dimensional description of scale interactions. Particular attention is devoted to the inner canopy, to identify the scales and structures responsible for sustaining the fluctuations, and to understand the mechanisms governing energy transfer from the shear layer. We also examine how canopy flexibility modulates these structures and their scale organisation.
The remainder of this paper is organised as follows. Section 2 outlines the physical set-up and the DNS database, and provides a brief overview of the AGKE framework. In §§ 3.1 and 3.2, we focus on the case with rigid canopies to identify characteristic flow scales and examine the mechanisms sustaining velocity fluctuations in the combined physical and scale space. The influence of canopy flexibility is then explored in § 3.3. Finally, concluding remarks and future perspectives are presented in § 4.
2. Mathematical formulation
2.1. Problem set-up and DNS database
We use the DNS database of turbulent flow over dense, submerged flexible and rigid canopies generated by Monti et al. (Reference Monti, Olivieri and Rosti2023). Full details of the numerical method and computational procedures are provided in Monti et al. (Reference Monti, Olivieri and Rosti2023) and Foggi Rota et al. (Reference Foggi Rota, Monti, Olivieri and Rosti2024), where the interested reader can also find validations against experimental measurements.
We study the turbulent flow in a planar half-channel containing a finite number of slender filaments clamped to the bottom wall (see figure 1). A Cartesian coordinate system is used, with
$x$
(
$u$
),
$y$
(
$v$
) and
$z$
(
$w$
) denoting the streamwise, wall-normal and spanwise directions (and corresponding velocity components). For simplicity, the alternative notation
$x_1$
(
$u_1$
),
$x_2$
(
$u_2$
) and
$x_3$
(
$u_3$
) is also occasionally adopted.
Sketch of the computational domain and of the reference frame.
The fluid motion is governed by the incompressible Navier–Stokes equations for a Newtonian fluid:
\begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla } \boldsymbol{\cdot }(\boldsymbol{u} \boldsymbol{u}) = - \frac {1}{\rho _f} \boldsymbol{\nabla }\! p + \nu {\nabla} ^2 \boldsymbol{u} + \boldsymbol{f}, \quad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} = 0, \end{equation}
where
$\boldsymbol{u}$
is the velocity vector,
$p$
is the pressure and
$\rho _f$
and
$\nu$
denote the fluid density and kinematic viscosity, respectively. The term
$\boldsymbol{f}$
represents the interaction force between the fluid and the filamentous structures. The bulk Reynolds number is fixed at
$\textit{Re}_b = U_b H / \nu = 5000$
for all cases, where
$U_b = (1/\varOmega ) \int _\varOmega u \, \mathrm{d} \varOmega$
is the bulk velocity,
$\varOmega$
is the volume of the computational domain and
$H$
is the channel half-height, measured from the wall to the centreline (see figure 1). The simulations maintain a constant
$U_b$
by suitably adapting in time the pressure gradient. The flexibility of the filaments is characterised by the dimensionless Cauchy number:
\begin{equation} \textit{Ca} = \frac {1}{2} \frac {\rho _f {\rm d} h^3 U_b^2}{\gamma }, \end{equation}
where
$d$
and
$h$
are the filament diameter and length, respectively, and
$\gamma$
is the bending stiffness. The Cauchy number quantifies the ratio between the fluid-induced deforming force and the elastic restoring force of the filaments, with larger values corresponding to more flexible structures.
The simulations were performed with the solver Fujin (Rosti Reference Rosti2026), which uses second-order finite-difference schemes in space and a second-order Adams–Bashforth scheme in time. A fractional-step method is employed for the temporal advancement of the solution. The computational domain extends for
$L_x = 2 \pi H$
,
$L_y = H$
and
$L_z = 1.5 \pi H$
. The adequacies of these values are discussed in Monti et al. (Reference Monti, Olivieri and Rosti2023) and further confirmations are provided by Sathe & Giometto (Reference Sathe and Giometto2024). No-slip and no-penetration boundary conditions are applied at the bottom wall (
$y = 0$
), symmetry conditions are used at the top boundary (
$y = H$
) and periodic conditions are enforced in the streamwise and spanwise directions. The computational domain is discretised with
$(N_x,N_y,N_z) = (1152,384,864)$
points. In the streamwise and spanwise directions a uniform distribution of points is adopted. In the wall-normal direction, within the canopy the distribution of points is uniform, while above the canopy tip the grid is linearly stretched up to the top boundary.
The canopy is formed by
$n_x \times n_z = 144 \times 108$
filaments in the streamwise and spanwise directions. Their length is
$h = 0.25H$
. To avoid preferential flow channelling, the surface of the solid wall is divided in a grid of
$n_x \times n_z$
equal squares; each is assigned one filament, whose location within the square is selected randomly. The fluid–solid interaction is dealt with using a Lagrangian immersed-boundary method (Yu Reference Yu2005; Huang et al. Reference Huang, Shin and Sung2007); see also Marchioli et al. (Reference Marchioli, Rosti and Verhille2025). Each filament is discretised with
$N_l$
Lagrangian points uniformly distributed along the stem length, with
$N_l=126$
for the rigid filaments (
$\textit{Ca} = 0$
) and
$N_l=32$
for the flexible ones (
$\textit{Ca} \neq 0$
). In the rigid case, a greater number of Lagrangian points is employed to resolve the larger velocity difference between the fluid and the structure. This denser distribution results in more uniform forcing, with roughly one Lagrangian point per Eulerian grid point. The no-slip condition
$\partial \boldsymbol{X}/\partial t = \boldsymbol{u}(\boldsymbol{X}(s,t),t)$
is enforced at each Lagrangian point (see the Appendix A). An anelastic collision kernel prevents geometrical overlap; however, within the present parameter range, contacts are extremely rare, and disabling the collision model produces statistically indistinguishable results. At every time step the fluid velocity
$\boldsymbol{u}$
is interpolated on the Lagrangian points, and the fluid–structure force acting on the filaments is evaluated as
\begin{equation} \boldsymbol{f}_{\!l} = - \beta \left ( \boldsymbol{u}_{l} - \frac {\partial \boldsymbol{X}}{\partial t} \right )\!, \end{equation}
where
$\partial \boldsymbol{X}/\partial t$
is the velocity of the filament,
$\boldsymbol{u}_{l}$
is the fluid velocity interpolated on the position of the filament and
$\beta$
is a properly chosen proportionality coefficient which is set equal to
$10$
. Note that we use the subscript ‘
$l$
’ for quantities defined in the Lagrangian frame of reference. The fluid–structure force
$\boldsymbol{f}_{\!l}$
is transferred to the Eulerian grid via a spreading function, producing the volume force
$\boldsymbol{f}$
appearing in (2.1). This approach enforces the no-slip and no-penetration conditions over a finite region rather than at isolated points. When the Lagrangian spacing is comparable with the Eulerian grid, overlapping kernels ensure continuous and effective enforcement along the filament.
The structural dynamics of the filaments for
$\textit{Ca} \neq 0$
is described with the generalised Euler–Bernoulli beam model:
\begin{equation} \begin{cases} \displaystyle \Delta \tilde {\rho } \frac {\partial ^2 \boldsymbol{X}}{\partial t^2} = \frac {\partial }{\partial s} \left ( T \frac {\partial \boldsymbol{X}}{\partial s} \right ) - \gamma \frac {\partial ^4 \boldsymbol{X}}{\partial s^4} - \boldsymbol{f}_{\!l}, \\[15pt] \displaystyle \frac {\partial \boldsymbol{X}}{\partial s} \boldsymbol{\cdot }\frac {\partial \boldsymbol{X}}{\partial s} = 1, \end{cases} \end{equation}
where
$\boldsymbol{X}(s,t)$
is the position of a point along the filament, a function of the curvilinear abscissa
$s$
and of the time
$t$
,
$\Delta \tilde {\rho }$
is the linear density difference between the filament and the fluid, defined as
$\Delta \tilde {\rho } = ( \rho _s - \rho _f )\pi d^2/4$
(
$\rho _s$
is the density per unit volume of the filaments), and
$T$
is the tension enforcing the inextensibility constraint. The lower end of the filaments (
$s = 0$
) is clamped to the bottom wall, while the other end (
$s = h$
) is free to oscillate, thus leading to the following set of boundary conditions:
$\boldsymbol{X}|_{s=0} = \boldsymbol{X}_0$
,
$\partial \boldsymbol{X}/\partial s|_{s=0} = (0,\phi _0,0)$
,
$\partial ^2\boldsymbol{X}/\partial s^2|_{s=h}=\boldsymbol{0}$
,
$\partial ^3\boldsymbol{X}/\partial s^3|_{s=h}=\boldsymbol{0}$
and
$T|_{s=h}=0$
, where
$\phi _0=0$
is the clamp angle of the stems. Equation (2.4) is solved following Huang et al. (Reference Huang, Shin and Sung2007), except for the bending term
$\gamma \partial ^4 \boldsymbol{X}/\partial s^4$
which is treated implicitly as in Banaei et al. (Reference Banaei, Rosti and Brandt2020) to allow for a larger time step.
In this study, we investigate three values of the Cauchy number, i.e.
$\textit{Ca} = 0$
,
$\textit{Ca} = 25$
and
$\textit{Ca} = 50$
, representing the rigid, intermediate and flexible regimes, respectively (see Monti et al. Reference Monti, Olivieri and Rosti2023). Figure 2 presents vertical profiles of the mean velocity, turbulent kinetic energy and Reynolds stresses for these cases.
For comparison, reference data from a planar open-channel flow without a canopy at the same bulk Reynolds number are also included in parts of the analysis (see § 3.1). The main computational parameters for all simulations are summarised in table 1.
Flow and canopy parameters for the rigid and flexible cases. Here
$\textit{Re}_b$
is the bulk Reynolds number,
$L_x$
and
$L_z$
are the streamwise and spanwise sizes of the computational domain,
$N_x$
,
$N_y$
and
$N_z$
are the numbers of grid points in the streamwise, wall-normal and spanwise directions,
$n_x$
and
$n_z$
are the numbers of stems in the streamwise and spanwise direction,
$\lambda$
is the canopy solidity,
$ \langle {\text{d}p/\text{d}x} \rangle$
is the average streamwise pressure gradient required to sustain the flow,
$ \langle {y_{\textit{tip}}} \rangle _f$
is the average wall-normal location of the tip of the filaments, with
$ \langle {\boldsymbol{\cdot }} \rangle _f$
denoting average in time and across filaments.
Table 1 Long description
A table comparing flow and canopy parameters for rigid and flexible cases. The table has 4 rows and 12 columns. Column headers are Ca, Re_b, (Lx, Lz)/H, N_x x N_y x N_z, n_x x n_z, λ, h/H, (dp/dx), (y_tip)_f, y_v0. Row 1: Rigid, 0, 5000, 2π, 1.5π, 1152 x 384 x 864, 144 x 108, 1.43, 0.25, 0.0441, 0.25, 186. Row 2: Intermediate, 25, 5000, 2π, 1.5π, 1152 x 384 x 864, 144 x 108, 1.43, 0.25, 0.0277, 0.19, 109. Row 3: Flexible, 50, 5000, 2π, 1.5π, 1152 x 384 x 864, 144 x 108, 1.43, 0.25, 0.0208, 0.135, 71. Row 4: Channel, –, 5000, 2π, 1.5π, 1152 x 384 x 864, –, –, –, 0.0034, –, –.
(a) Mean streamwise velocity
$U$
, (b) turbulent kinetic energy
$k$
and (c–e) Reynolds stresses
$ \langle {uu} \rangle$
,
$ \langle {vv} \rangle$
and
$ \langle {ww} \rangle$
for
$\textit{Ca} = 0$
, 25 and 50. All quantities are non-dimensionalised by
$U_b$
and
$H$
. The dashed lines indicate the average wall-normal location of the tip of the filaments, i.e.
$ \langle {y_{\textit{tip}}} \rangle _f/H = 0.25, 0.19, 0.135$
for
$\textit{Ca}=0, 25, 50$
, respectively;
$ \langle {\boldsymbol{\cdot }} \rangle _f$
denotes average in time and across filaments. The same colour code is used for solid and dashed lines.
Figure 2 Long description
Panel A: A line graph shows the mean streamwise velocity (U/U_b) against the normalized wall-normal position (y/H). The graph includes three lines representing different capillary numbers (Ca = 0, Ca = 25, Ca = 50) in blue, red, and yellow respectively. An inset zooms in on the velocity near the top of the graph. Panel B: A line graph displays turbulent kinetic energy (k/U_b^2) against the normalized wall-normal position (y/H) for the same capillary numbers. Panel C: A line graph illustrates the Reynolds stress component (uu/U_b^2) against the normalized wall-normal position (y/H) for the same capillary numbers. Panel D: A line graph shows the Reynolds stress component (vv/U_b^2) against the normalized wall-normal position (y/H) for the same capillary numbers. Panel E: A line graph depicts the Reynolds stress component (ww/U_b^2) against the normalized wall-normal position (y/H) for the same capillary numbers. Dashed lines indicate the average wall-normal location of the tip of the filaments for each capillary number.
For all the analysed cases, a total of 50 flow-field snapshots, distributed over a time of
$30 H/U_b$
, are employed for computing the flow statistics. After the Reynolds decomposition, the velocity field is written as the sum of the mean flow
$\boldsymbol{U}= \left \langle {\boldsymbol{u}} \right \rangle$
and the fluctuations, where the
$ \left \langle {\boldsymbol{\cdot }} \right \rangle$
operator represents the statistical average along the homogeneous directions and in time. For simplicity, hereinafter we use uppercase letters for mean quantities and lowercase letters for fluctuations.
Unless stated otherwise, all quantities are non-dimensionalised using the bulk velocity
$U_b$
and the channel half-height
$H$
. Outer units are preferred for our analysis, as they directly relate to the large-scale motions in the unobstructed region of the channel (Peruzzi et al. Reference Peruzzi, Poggi, Ridolfi and Manes2020). We also use the superscript ‘
$+$
’ to denote normalisation in viscous (inner) units, where quantities are made dimensionless using the friction velocity
$u_\tau$
and the kinematic viscosity
$\nu$
. In canopy flows, when viscous units are applied to quantities in the outer layer, the friction velocity is evaluated at the virtual origin
$y_{vo}$
, defined as the effective origin of the logarithmic velocity profile. Following Foggi Rota et al. (Reference Foggi Rota, Monti, Olivieri and Rosti2024), this virtual origin is determined by matching the mean velocity profile above the canopy to a canonical logarithmic law, with the corresponding friction velocity computed as
$u_\tau = ( \nu \mathrm{d}U/\mathrm{d}y - \langle u v \rangle - \int ^H_y \left \langle {f_x} \right \rangle (s) \,\mathrm{d}s )^{1/2}|_{y=y_{vo}}$
. For the smooth-wall reference case, the friction velocity is instead evaluated directly at the physical wall (
$y=0$
). Notably, the estimated values of
$u_\tau$
and
$y_{vo}$
are highly robust; applying alternative identification strategies, such as the diagnostic-function matching proposed by Chen & García-Mayoral (Reference Chen and García-Mayoral2023), yields only marginal differences that do not alter any physical results. Finally, to eliminate any potential coordinate ambiguity, all wall-normal positions throughout this paper are consistently expressed as distances from the physical wall (
$y$
), regardless of whether they are normalised in inner or outer units. For reference, we report the evaluated values of
$y^+_{vo}$
for the different configurations in table 1.
2.2. The anisotropic generalised Kolmogorov equations
The AGKE, introduced by Gatti et al. (Reference Gatti, Chiarini, Cimarelli and Quadrio2020) and later extended by Gattere et al. (Reference Gattere, Chiarini, Gallorini and Quadrio2023), provide the exact budget equations for the components of the second-order velocity structure function tensor,
$\langle \delta u_i \delta u_{\!j} \rangle$
. Here,
$\delta u_i = u_i(\boldsymbol{x}_\beta ,t) - u_i(\boldsymbol{x}_\alpha ,t)$
denotes the increment of the fluctuating velocity component
$u_i$
between two points
$\boldsymbol{x}_\alpha$
and
$\boldsymbol{x}_\beta$
. The midpoint
$\boldsymbol{X} \equiv (\boldsymbol{x}_\alpha + \boldsymbol{x}_\beta )/2$
and the separation vector
$\boldsymbol{r} \equiv \boldsymbol{x}_\beta - \boldsymbol{x}_\alpha$
define the position and the scale, respectively. In the general case, the structure function tensor
$\langle \delta u_i \delta u_{\!j} \rangle$
depends on seven independent variables: the three components of
$\boldsymbol{X}$
, the three components of
$\boldsymbol{r}$
and time
$t$
.
The tensor
$\langle \delta u_i \delta u_{\!j} \rangle$
combines single-point and two-point statistics. Specifically, it can be expressed as
\begin{equation} \langle \delta u_i \delta u_{\!j} \rangle (\boldsymbol{X}, \boldsymbol{r}, t) = \mathsf{V}_{ij}(\boldsymbol{X}, \boldsymbol{r}, t) - \mathsf{R}_{ij}(\boldsymbol{X}, \boldsymbol{r}, t) - \mathsf{R}_{ij}(\boldsymbol{X}, -\boldsymbol{r}, t), \end{equation}
where
$\mathsf{R}_{ij}(\boldsymbol{X}, \boldsymbol{r}, t) = \langle u_i(\boldsymbol{x}_\alpha ,t) u_{\!j}(\boldsymbol{x}_\beta ,t) \rangle$
is the two-point velocity correlation tensor and
$\mathsf{V}_{ij} = \langle u_i u_{\!j} \rangle |_{\boldsymbol{x}_\alpha } + \langle u_i u_{\!j} \rangle |_{\boldsymbol{x}_\beta }$
is the sum of the single-point Reynolds stresses at the two
$\boldsymbol{x}_\alpha$
and
$\boldsymbol{x}_\beta$
points. The structure function thus reflects the combined effect of local turbulence intensity and the correlation of velocity fluctuations across scales. Notably,
$\mathsf{R}_{ij}$
carries the essential scale-dependent information.
A scalar quantity of interest is
$\langle \delta q^2 \rangle = \sum _i \langle \delta u_i \delta u_i \rangle$
, often interpreted as a proxy for the turbulent kinetic energy associated with eddies of size
$r = |\boldsymbol{r}|$
at position
$\boldsymbol{X}$
, i.e. the scale energy. However, this interpretation is only approximate:
$\langle \delta q^2 \rangle$
reflects the cumulative energy of all eddies with scales up to
$r$
, rather than the energy at scale
$r$
itself. In the idealised case of homogeneous isotropic turbulence, where statistics depend only on
$r$
, this becomes evident from the identity
$\langle \delta q^2 \rangle = 2 \mathsf{V}_{ii} - 2 \mathsf{R}_{ii}(r)$
, so that
$\langle \delta q^2 \rangle = 0$
at
$r=0$
and tends to
$2 \mathsf{V}_{ii}$
as
$r \to \infty$
. However, as discussed by Davidson & Pearson (Reference Davidson and Pearson2005), this behaviour also includes contributions from large-scale eddies due to their associated enstrophy, and thus cannot be regarded as purely local in scale.
The derivation of the AGKE begins with the equation for the velocity difference
$\delta u_i$
, obtained by subtracting the Navier–Stokes equations at points
$\boldsymbol{x}_\alpha$
and
$\boldsymbol{x}_\beta$
. The procedure involves multiplication by
$\delta u_{\!j}$
, expressing all terms as two-point differences or sums, transforming the coordinate system from
$(\boldsymbol{x}_\alpha , \boldsymbol{x}_\beta )$
to
$(\boldsymbol{X}, \boldsymbol{r})$
and averaging. Full details are provided in Gatti et al. (Reference Gatti, Chiarini, Cimarelli and Quadrio2020) and Gattere et al. (Reference Gattere, Chiarini, Gallorini and Quadrio2023). The AGKE describe the scale-by-scale balance of Reynolds stresses, linking the temporal evolution of
$\langle \delta u_i \delta u_{\!j} \rangle$
at a given scale and position to the combined effects of production, transport and dissipation. In the present study, the flow is statistically stationary and homogeneous in two directions, reducing the number of independent variables to four:
$(Y, r_x, r_y, r_z)$
. The AGKE tailored to this configuration read
\begin{equation} \begin{gathered} \frac {\partial }{\partial r_k} \biggl ( \underbrace { \langle \delta u_i \delta u_{\!j} \delta U_k \rangle }_{\phi ^{m}_{k,ij}} + \underbrace { \langle \delta u_i \delta u_{\!j} \delta u_k \rangle }_{\phi ^{t}_{k,ij}} - \underbrace {2\nu \frac {\partial }{\partial r_k} \langle \delta u_i \delta u_{\!j} \rangle }_{\phi ^{v}_{k,ij}} \biggr ) \\+ \frac {\partial }{\partial Y} \biggl ( \underbrace { \langle v^* \delta u_i \delta u_{\!j} \rangle }_{\psi ^{t}_{ij}} + \underbrace { \frac {1}{\rho _f} \langle \delta p \delta u_{\!j} \rangle \delta _{i2} + \frac {1}{\rho _f} \langle \delta p \delta u_i \rangle \delta _{j2} }_{\psi ^{p}_{ij}} - \underbrace { \frac {\nu }{2} \frac {\partial }{\partial Y} \langle \delta u_i \delta u_{\!j} \rangle }_{\psi ^{v}_{ij}} \biggr ) \\= \underbrace {- \langle v^* \delta u_{\!j} \rangle \delta \left ( \frac {\text{d}U}{\text{d}y} \right ) \delta _{i1} - \langle v^* \delta u_i \rangle \delta \left ( \frac {\text{d}U}{\text{d}y} \right ) \delta _{j1} - \langle \delta v \delta u_{\!j} \rangle \left ( \frac {\text{d}U}{\text{d}y} \right )^* \delta _{i1} - \langle \delta v \delta u_i \rangle \left ( \frac {\text{d}U}{\text{d}y} \right )^* \delta _{j1}}_{ {P} _{ij}} \\+ \underbrace { \frac {1}{\rho _f} \biggl \langle \delta p \frac {\partial \delta u_i}{\delta X_{\!j}} \biggr \rangle + \frac {1}{\rho _f} \biggl \langle \delta p \frac {\partial \delta u_{\!j}}{\delta X_i} \biggr \rangle }_{\varPi _{ij}} - \underbrace {4\epsilon ^*_{ij}}_{ {D} _{ij}} + \underbrace { \langle \delta u_{\!j} \delta f_i \rangle + \langle \delta u_i \delta f_{\!j} \rangle }_{ {F} _{ij}}, \end{gathered} \end{equation}
where
$i = 1,2,3$
and
$j=1,2,3$
. Here
$\delta _{ij}$
is the Kronecker delta, the asterisk denotes the average of a given quantity between the two points
$\boldsymbol{x}_\alpha$
and
$\boldsymbol{x}_\beta$
and
$\epsilon _{ij} = \nu \langle {\partial u_i/\partial x_k \ \partial u_{\!j}/\partial x_k} \rangle$
is the pseudo-dissipation tensor. The
$\boldsymbol{\phi }^{m}_{ij}$
,
$\boldsymbol{\phi }^{t}_{ij}$
and
$\boldsymbol{\phi }^{v}_{ij}$
terms are the mean, turbulent and viscous fluxes in the space of scales. The
${\psi }^{t}_{ij}$
,
${\psi }^{p}_{ij}$
and
${\psi }^{v}_{ij}$
terms, instead, are the turbulent, pressure and viscous fluxes in the physical space; they describe how energy is transported in the wall-normal direction. In addition,
$ {P} _{ij}$
,
$\varPi _{ij}$
,
$ {D} _{ij}$
and
$ {F} _{ij}$
are the production, pressure–strain, dissipation and forcing terms, respectively. All the listed terms are named in analogy to the single-point budgets (see e.g. Pope (Reference Pope2000). The production term
$ {P} _{ij}$
couples the mean and fluctuating fields and describes the scales and positions at which they exchange energy; the pressure–strain term
$\varPi _{ij}$
couples pressure and velocity derivatives, and describes the energy redistribution among the three normal stresses; the forcing term
$ {F} _{ij}$
describes the effect of the interaction between fluid and filaments in the combined space of scales and positions. Summing together the terms on the right-hand side yields the source term
$\xi _{ij} = {P} _{ij} + \varPi _{ij} + {D} _{ij} + {F} _{ij}$
. The sum of the two-point budgets for the three normal stresses leads to the generalised Kolmogorov equation (Hill Reference Hill2002).
In this work, the terms of the AGKE are computed using a computer code, written in the CPL language (Luchini Reference Luchini2021), and originally developed for the isotropic version of the AGKE by Gatti et al. (Reference Gatti, Remigi, Chiarini, Cimarelli and Quadrio2019). For maximum computational efficiency, whenever possible, the code computes two-point correlations in the spectral domain, leveraging Parseval’s theorem.
3. Results
We investigate the mechanisms sustaining turbulence in dense canopy flows by analysing the dynamics both in physical space and across scales. We first identify in § 3.1 the dominant flow scales and their associated coherent structures. Building on this framework, § 3.2 examines the production, pressure–strain redistribution and inter-scale transfer of velocity fluctuations, assessing how these mechanisms vary in space and how they relate to the underlying flow structures. The rigid canopy configuration (
$\textit{Ca} = 0$
) is adopted as the reference case, while the effect of filament flexibility is discussed separately in § 3.3. For comparison, data from a smooth open-channel flow at matched
$\textit{Re}_b$
are included in § 3.1 to isolate and highlight the impact of the canopy on the dominant scales and coherent motions.
3.1. Flow structure
In this section, we identify the dominant scales in the flow and relate them to coherent structures. Structure functions and two-point correlations are first employed to quantify the characteristic velocity scales. Conditional averaging and filtering are then used to associate these scales with the underlying motions, including KH-like rollers and hairpin vortices. This approach combines objective, scale-resolved diagnostics with physically interpretable representations of the coherent dynamics.
3.1.1. Structure functions and flow scales
The canopy flow exhibits a broad range of scales, with turbulence largely confined to the outer layer and the interior remaining relatively quiescent (figure 3). Peak turbulent activity occurs near the canopy tips, where fluid–structure interactions generate intense shear layers. Coherent structures observed include low- (LSS) and high-speed streaks (HSS) of streamwise velocity, along with large-scale vortical motions aligned in both the streamwise and spanwise directions (see also Finnigan Reference Finnigan2000).
Coherent structures in a snapshot of the canopy flow (
$\textit{Ca}=0$
). Isosurfaces of streamwise velocity fluctuations at
$u = \pm 0.6$
are coloured in red and blue; isosurfaces of negative filtered pressure fluctuations
$\tilde {p} = -0.06$
(see § 3.1.2) are in cyan. The background planes plot streamwise velocity fluctuations
$u$
in greyscale.
Figure 4 quantifies the dominant scales through the diagonal components of the structure function tensor
$\langle \delta u_i \delta u_{\!j} \rangle$
in the
$r_y = 0$
plane, with smooth-wall open-channel data at matched
$\textit{Re}_b$
provided for reference. Peaks in these distributions identify characteristic length scales and highlight the wall-normal variation of the fluctuations.
Structure functions
$\langle \delta u_i \delta u_{\!j} \rangle$
evaluated in the
$r_y = 0$
space. (a,c,e) Canopy flow (
$\textit{Ca}=0$
); (b,d, f) smooth open-channel flow. From top to bottom, rows correspond to
$\langle \delta u \delta u \rangle$
,
$\langle \delta v \delta v \rangle$
and
$\langle \delta w \delta w \rangle$
, respectively. Grey isosurfaces highlight regions of intense structure function magnitude, corresponding to local peaks of
$\langle \delta u_i \delta u_{\!j} \rangle$
observed in the
$r_y = r_z = 0$
and
$r_x = r_y = 0$
planes, and associated with negative values of
$\mathsf{R}_{ij}$
at those scales. The values are reported in the figures. The absence of emerging isosurfaces from these planes indicates that no pronounced local peaks of
$\langle \delta u_i \delta u_{\!j} \rangle$
and
$\mathsf{R}_{ij}$
are present at the corresponding scales.
Figure 4 Long description
Panel A: A 3D contour plot shows the structure function for canopy flow in the uu plane. The x-axis is labeled rz, the y-axis is labeled Y, and the z-axis is labeled rx. The color scale ranges from 0 to 0.2. Grey isosurfaces highlight regions of intense structure function magnitude. Panel B: A 3D contour plot shows the structure function for smooth open-channel flow in the uu plane. The x-axis is labeled rz, the y-axis is labeled Y, and the z-axis is labeled rx. The color scale ranges from 0 to 0.05. Grey isosurfaces highlight regions of intense structure function magnitude. Panel C: A 3D contour plot shows the structure function for canopy flow in the vv plane. The x-axis is labeled rz, the y-axis is labeled Y, and the z-axis is labeled rx. The color scale ranges from 0 to 0.08. Grey isosurfaces highlight regions of intense structure function magnitude. Panel D: A 3D contour plot shows the structure function for smooth open-channel flow in the vv plane. The x-axis is labeled rz, the y-axis is labeled Y, and the z-axis is labeled rx. The color scale ranges from 0 to 0.006. Grey isosurfaces highlight regions of intense structure function magnitude. Panel E: A 3D contour plot shows the structure function for canopy flow in the ww plane. The x-axis is labeled rz, the y-axis is labeled Y, and the z-axis is labeled rx. The color scale ranges from 0 to 0.12. Grey isosurfaces highlight regions of intense structure function magnitude. Panel F: A 3D contour plot shows the structure function for smooth open-channel flow in the ww plane. The x-axis is labeled rz, the y-axis is labeled Y, and the z-axis is labeled rx. The color scale ranges from 0 to 0.01. Grey isosurfaces highlight regions of intense structure function magnitude.
Turbulent activity in the outer layer peaks near the canopy tips, where the diagonal components of the structure function tensor reach their maximum, and gradually decays towards the channel centre. At this Reynolds number, turbulence in the core is weak, and remains strongly attenuated within the canopy. Comparison with smooth-wall channel flow reveals both similarities and differences: the streamwise and wall-normal structure functions exhibit distinct peaks at
$(Y, r_x, r_z) = (0.273, 0, 0.327)$
and
$(0.365, 0, 0.284)$
for the canopy, and at
$(0.049, 0, 0.218)$
and
$(0.209, 0, 0.262)$
for the smooth-wall channel, corresponding in viscous units to
$(Y^+, r_x^+, r_z^+) = (257, 0, 308)$
and
$(344, 0, 268)$
for the canopy, and
$(14, 0, 64)$
and
$(61, 0, 77)$
for the smooth wall. These values indicate coherent LSS and HSS and streamwise-aligned vortices in both flows, consistent with observations in high-
$Re$
flexible canopies (Löhrer & Fröhlich Reference Löhrer and Fröhlich2025). Peaks in the structure functions correspond to minima in the two-point correlation functions
$\mathsf{R}_{ij}$
, with negative
$\mathsf{R}_{22}$
at
$r_z \ne 0$
reflecting upwash and downwash motions induced by streamwise-aligned vortices. In canopy flows, these motions are associated with the legs of hairpin vortices, which populate the region just above the tip of the filaments; see Finnigan et al. (Reference Finnigan, Shaw and Patton2009) and § 3.1.2. The coherence of these fluctuations is quantified via correlation coefficients
$\rho _{ij}(Y,r_x,r_z) = \mathsf{R}_{ij}(Y,r_x,r_z)/ \langle { u_i u_{\!j} } \rangle (Y)$
at the structure function peaks, yielding
$\rho _{11}=-0.168$
and
$\rho _{22}=-0.04$
for the canopy, and
$\rho _{11}=-0.150$
and
$\rho _{22}=-0.131$
for the smooth wall. This suggests that LSS/HSS are more coherent in canopy flows than in canonical channel flows, while the opposite is true for streamwise-aligned vortices.
A clear distinction between the outer-layer canopy flow and the smooth-wall channel is evident in the
$r_z = 0$
plane (figure 4
c,d). In the canopy case, the map of
$\langle \delta v \delta v \rangle$
displays a secondary, weaker peak at
$Y \approx 0.365$
(
$Y^+ \approx 344$
) and streamwise separation
$r_x = 0.436$
(
$r^+_x = 411$
), which is absent in the smooth-wall configuration. This feature corresponds to spanwise-aligned KH-type vortices originating from the shear-layer instability near the canopy tip (see § 3.1.2), which induce vertical motions along their sides, producing negative
$\mathsf{R}_{22}$
for
$r_x \ne 0$
. The associated correlation coefficient,
$\rho _{22}=-0.024$
, indicates that these KH structures are less coherent and statistically weaker than the primary LSS/HSS and streamwise-aligned vortices.
Contours of
$ \langle {\delta u \delta u} \rangle$
(a) and
$\mathsf{R}_{11}$
(b) in the
$(Y,r_z)$
space for
$\textit{Ca}=0$
. The white crosses indicate the location of the maxima. The white crosses with a circle indicate the location of the secondary maxima. In (b) the black solid line represents the
$\mathsf{R}_{11}=0$
contour. The blue dashed line is the wall-normal position of the canopy tip. (c) Schematic representation of the process which leads to the formation of two couples of streaks.
A closer inspection of
$ \left \langle {\delta u \delta u} \right \rangle$
in the
$r_x = r_y = 0$
plane (figure 5) reveals a secondary peak at larger spanwise separations (
$r_z \approx 1.6$
,
$r_z^+ \approx 1508$
) in the canopy, absent in the smooth-wall case. At
$Y = 0.28$
(
$Y^+ = 264$
), the corresponding correlation coefficient is
$\rho _{11} = -0.018$
. This secondary peak reflects paired HSS and LSS, resulting from a spanwise modulation of KH rolls induced by the mean shear near the canopy tips, as illustrated schematically in figure 5(c) and discussed in Finnigan (Reference Finnigan2000), Finnigan et al. (Reference Finnigan, Shaw and Patton2009), Bailey & Stoll (Reference Bailey and Stoll2016) and Tschisgale et al. (Reference Tschisgale, Löhrer, Meller and Fröhlich2021). Peaks in
$ \left \langle {\delta w \delta w} \right \rangle$
in the
$r_x = 0$
plane further confirm the lateral organisation associated with these KH-induced streaks. This mechanism contrasts with smooth-wall turbulence, where streaks are primarily driven by vertical motions from quasi-streamwise vortices (Jiménez Reference Jiménez2018), and no secondary peak in
$ \left \langle {\delta u \delta u} \right \rangle$
or
$\mathsf{R}_{11}$
is observed.
Organisation of velocity fluctuations within the canopy layer (
$\textit{Ca}=0$
). (a–c) Structure functions
$ \langle {\delta u_i \delta u_{\!j}} \rangle$
shown in the
$r_y=0$
space, zoomed inside the canopy region
$Y\lt 0.2$
. (d) Profiles of
$ \langle { \delta u_i \delta u_{\!j} } \rangle$
along the blue dashed lines indicated in (a,c).
Within the inner canopy layer (
$Y \leqslant 0.2$
; figure 6), velocity fluctuations are dominated by the streamwise and spanwise components, consistent with a quasi-two-dimensional flow imposed by the impermeable wall. As expected,
$u$
fluctuations are stronger than
$w$
fluctuations. Despite overall turbulence suppression, weak but discernible peaks remain in the wall-parallel structure functions, reflecting residual scale organisation:
$ \left \langle {\delta u \delta u} \right \rangle$
peaks at
$r_x \approx 0.4$
–
$0.6$
,
$r_z = 0$
, while
$ \left \langle {\delta w \delta w} \right \rangle$
peaks at
$r_x = 0$
,
$r_z \approx 0.3$
–
$0.5$
. The corresponding correlation coefficients,
$\rho _{11} = \rho _{33} = -0.2$
, indicate that these structures are statistically significant and dominate the local fluctuations relative to the quiescent background. This organisation reflects the modulation imposed by scale-dependent motions above the canopy, as detailed in § 3.2. Streamwise- and spanwise-aligned vortices generate intense sweep events that induce localised downwash into the canopy. These vertical motions are redirected into wall-parallel components via pressure–strain interactions at the wall: sweeps associated with streamwise-aligned vortices enhance spanwise fluctuations, while those associated with spanwise-aligned vortices enhance streamwise fluctuations, reproducing the characteristic scales observed in
$ \left \langle {\delta w \delta w} \right \rangle$
and
$ \left \langle {\delta u \delta u} \right \rangle$
, respectively.
(a) Streamwise and (b) spanwise scales at which
$\langle \delta u_i \, \delta u_{\!j} \rangle$
attain their maximum values as a function of the wall-normal location
$Y$
for
$\textit{Ca}=0$
. The dashed blue line indicates the position of the canopy tip. The symbol size corresponds to the correlation coefficient, with larger (smaller) symbols indicating higher (lower) correlation values.
Figure 7 summarises the characteristic length scales as a function of
$Y$
, along with the corresponding correlation coefficients. Near the canopy tips, structures are compact; moving into the canopy or towards the outer layer, scales increase, reflecting a broadening of influence with distance. Correlation coefficients peak near the wall, demonstrating that the flow structures above the canopy imprint a coherent organisation on the otherwise quiescent interior.
3.1.2. Eduction of the flow structures
In this section, we complement the statistically robust analysis of structure functions with the identification of the coherent flow structures. The goal is to relate the characteristic scales identified previously to specific flow features, providing a qualitative interpretation of their spatial organisation. While the previous section delivered quantitative information, the present analysis uses conditional averaging and filtering to isolate representative structures, yielding primarily qualitative insight.
The KH-like rolls are extracted from the instantaneous flow field by filtering out small-scale fluctuations using a Fourier filter, with filtered quantities denoted by a tilde. A Fourier transform is applied in the homogeneous streamwise and spanwise directions, retaining modes with wavenumbers
$\kappa _x \leqslant \kappa _{x,{th}} = 12$
and
$\kappa _z \leqslant \kappa _{z,{th}} = 3$
, chosen to capture the characteristic streamwise spacing and large-scale spanwise modulation of the rolls (Raupach et al. Reference Raupach, Finnigan and Brunei1996; Bailey & Stoll Reference Bailey and Stoll2016). Moderate variations in the cutoffs do not qualitatively alter the results. The filtered field retains approximately
$5\,\%$
of the wall-normal kinetic energy, providing a measure of the rolls’ energetic footprint. Figure 8 illustrates the filtered flow. Isosurfaces of negative pressure
$\tilde {p}$
(figure 8
a) and positive
$Q$
(figure 8
b) identify the KH roll cores, confirming spanwise-coherent structures above the canopy. Contours of
$\tilde {u}$
and
$\tilde {v}$
show a clear correlation between large-scale streamwise and vertical motions, linking the rolls to sweep and ejection events at the canopy tips and facilitating momentum exchange across the canopy interface.
A qualitative characterisation is obtained via conditional averaging. Events are defined as local minima in the filtered pressure,
$\tilde {p} \leqslant \tilde {p}_{{th}} = -0.065$
, within
$y \in (0.25, 0.4)$
; conditionally averaged quantities are denoted by
$\hat {\boldsymbol{\cdot }}$
. Small variations in
$\tilde {p}_{{th}}$
do not qualitatively affect the results. The resulting KH roll (figure 8
c) accounts for roughly
$3\,\%$
of the wall-normal kinetic energy. The negative
$\hat {p}$
contour confirms spanwise coherence, while the conditional velocity field exhibits upward motion at the rear and downward motion at the front, consistent with negative spanwise vorticity. Despite the dependence on the chosen
$\hat {p}$
contour, the structure shows a streamwise length
$\ell _x \approx 0.4{-}0.5$
(
$\ell _x^+ \approx 377{-}471$
), in agreement with the scale inferred from the
$ \left \langle {\delta v \delta v} \right \rangle$
map at
$r_y = r_z = 0$
.
Instantaneous and conditionally averaged KH rolls for
$\textit{Ca} = 0$
. (a) Isosurfaces of pressure fluctuations
$\tilde {p} = -0.04$
, coloured by the wall-normal velocity
$\tilde {v}$
. The horizontal slice shows contours of
$\tilde {u}\tilde {v} \lt 0$
at
$y \approx h$
. (b) Contours of the
$Q$
-criterion (
$\tilde {Q} = 1.5$
) computed from the complete (mean plus fluctuating) Fourier-filtered field, coloured by the spanwise vorticity
$\tilde {\omega }_z$
. In both panels, the mean flow is directed from left to right. (c) Conditionally averaged KH roll. The contour of negative pressure
$\hat {p} = -0.065$
is shown, coloured by the wall-normal velocity
$\hat {v}$
.
The educed hairpin-like vortices for
$\textit{Ca}=0$
. (a) Head-up hairpin-like eddies; (b) head-down hairpin-like eddies. The isosurface corresponds to
$\hat {Q} = 1$
, coloured with streamwise vorticity
$\hat {\omega }_x$
. The streamlines of conditionally averaged two-dimensional velocity, i.e.
$\boldsymbol{\hat {u}}_{2D} = (\hat {v}, \hat {w})$
, are coloured with the wall-normal velocity component
$\hat {v}$
.
Figure 9 Long description
Panel A: A 3D visualization showing head-up hairpin-like eddies. The isosurface is colored with streamwise vorticity. Streamlines of conditionally averaged two-dimensional velocity are colored with the wall-normal velocity component. Panel B: A 3D visualization showing head-down hairpin-like eddies. The isosurface is colored with streamwise vorticity. Streamlines of conditionally averaged two-dimensional velocity are colored with the wall-normal velocity component.
Hairpin-like vortices often emerge from the breakdown of KH rolls under mean shear (Finnigan et al. Reference Finnigan, Shaw and Patton2009), although they can also form via longitudinal stretching of weak vorticity between adjacent KH rolls (Corcos & Lin Reference Corcos and Lin1984; Neu Reference Neu1984). The resulting vortices feature streamwise-aligned, counter-rotating legs. We identify such events using a conditional averaging strategy analogous to that for KH rolls, applied to the unfiltered velocity field at wall-normal distances
$y \in (0.25, 0.4)$
. Hairpin-like eddies are detected as local minima of the shear stress,
$uv \lt uv_{th}$
(Bailey & Stoll Reference Bailey and Stoll2016), with the threshold set to
$uv_{th}/u^{*2} = -12$
, where
$u^{*2} = - \min _{y \in [0,H]} \left \langle {uv} \right \rangle$
. Variations in the threshold within
$uv_{th}/u^{*2} \in [-14, -10]$
produce qualitatively similar results. This method captures vortices that transport low- or high-momentum fluid between their legs depending on rotation, generating strong negative
$uv$
events (Brunet Reference Brunet2020; Shig et al. Reference Shig, Babin, Shnapp, Fattal, Liberzon and Bohbot-Raviv2023). Vortices are classified as head-up (
$v \gt 0$
, ejections) or head-down (
$v \lt 0$
, sweeps) based on the sign of
$v$
at the
$uv$
minimum. Figure 9 shows the resulting conditionally averaged hairpin-like structures. The geometry clearly displays the characteristic counter-rotating legs. The characteristic spanwise scale of their legs is
$\ell _z \approx 0.25$
(
$\ell _z^+ \approx 236$
), consistent with the dominant scale inferred from the
$ \left \langle {\delta v \delta v} \right \rangle$
map at
$r_x = r_y = 0$
. These conditional fields contain approximately
$2.5\,\%$
of the vertical kinetic energy. Some sensitivity of
$\ell _z$
to the selected
$Q$
-criterion contour level is noted.
To place the energetic footprint of the KH rolls and hairpin-like structures in context, we compare them with the quasi-streamwise vortices of a canonical turbulent plane channel flow at
$\textit{Re}_\tau = 180$
, isolated using the eduction procedure of Gallorini et al. (Reference Gallorini, Quadrio and Gatti2022). The conditionally averaged quasi-streamwise vortices account for approximately
$7\,\%$
of the total turbulent kinetic energy and
$5.5\,\%$
of its wall-normal component – values directly comparable to the global percentages obtained for the KH rolls and hairpin-like structures in the present canopy flow. These seemingly modest global fractions do not imply limited dynamical significance; indeed, two-point correlations and scale-by-scale energy analyses confirm that the identified structures in both flows actively mediate vertical momentum transport and inter-scale energy transfer. Furthermore, because wall-bounded turbulence is strongly inhomogeneous in the wall-normal direction, global energy normalisation inherently masks the intensity of localised coherent motions. When the energetic contribution of the KH rolls is re-evaluated locally – using spatial integrals restricted to horizontal planes or thin three-dimensional slabs (
$\Delta y \leqslant 0.1H$
) centred around the vortex cores – their contribution to the local wall-normal kinetic energy increases to approximately
$8\,\%$
.
3.2. Sustaining mechanisms
Having identified the dominant scales, we now examine the production, redistribution and transfer of turbulent kinetic energy towards dissipative scales, together with the energetic coupling between the outer flow and the canopy. We also analyse the mechanisms shaping the structure of the velocity fluctuations within the canopy.
To facilitate the analysis, we focus on the reduced space at
$r_y = 0$
, which retains the essential physics while simplifying the problem. In this space, the governing equations preserve their conservative form:
\begin{equation} \frac {\partial \phi _{x,ij}}{\partial r_x} + \frac {\partial \phi _{z,ij}}{\partial r_z} + \frac {\partial \psi _{ij}}{\partial Y} = \zeta _{ij}(r_x, r_z, Y), \end{equation}
demonstrating that in the
$r_y=0$
space the evolution of the fluxes is driven by the modified source term:
\begin{equation} \begin{aligned} \zeta _{ij} = & - \langle \delta v \delta u_{\!j} \rangle \left ( \frac {\mathrm{d}U}{\mathrm{d}y} \right ) \delta _{i1} - \langle \delta v \delta u_i \rangle \left ( \frac {\mathrm{d}U}{\mathrm{d}y} \right ) \delta _{j1} + \frac {1}{\rho _f} \left\langle \delta p \frac {\partial \delta u_i}{\partial X_{\!j}} \right\rangle + \frac {1}{\rho _f} \left\langle \delta p \frac {\partial \delta u_{\!j}}{\partial X_i} \right\rangle \\ & + \langle \delta u_i \delta f_{\!j} \rangle + \langle \delta u_{\!j} \delta f_i \rangle - 4 \epsilon _{ij} - \frac {\partial \phi _{y,ij}}{\partial r_y}. \end{aligned} \end{equation}
We analyse the scale energy budget,
$\langle \delta q^2 \rangle = \sum _i \langle \delta u_i \delta u_i \rangle$
, to identify the scales and locations of energy production and dissipation, and to highlight inter-scale transfers. The source term,
$\xi = \sum _i \xi _{ii}$
, captures the scale-dependent balance between production and dissipation, while the flux terms quantify transfers across scales and physical space. Near-wall pressure–strain redistribution and interfacial energy transfer are analysed to elucidate the mechanisms underlying the coherent motions observed within the canopy layer (§ 3.1).
3.2.1. Sources and sinks
(a) Colour map of the source term
$\xi$
in the
$r_y = 0$
plane for
$\textit{Ca} = 0$
. Field lines of the flux vector
$\boldsymbol{\varPhi } = (\phi _x, \phi _z, \psi )$
are shown in the same plane, coloured by the divergence
$\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\varPhi }$
. The grey surface represents the contour
$\xi = 0.9\xi _{max }$
. Evolution of the values of
$r_x$
,
$r_z$
and
$Y$
along a representative field line of
$\boldsymbol{\varPhi }$
for (b) group I and (c) group II. The dimensionless arch length
$s$
is defined as
$s=1/s_{max} \int ^{s_{max}}_{0} \text{d}s$
, with
$\text{d}s = \sqrt {\text{d}r^2_x + \text{d}r^2_z + \text{d}Y^2}$
. The red background is for
$\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\varPhi } \lt 0$
, and the blue one for
$\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\varPhi } \gt 0$
. Vectors of (d, f) turbulent transport
$\boldsymbol{\varPhi }^{t}$
and (e,g) viscous diffusion
$\boldsymbol{\varPhi }^{v}$
in the
$r_y = 0$
plane, shown at different wall-normal locations for the rigid canopy. Panels (d,e) correspond to
$Y = 0.15H$
(within the canopy), while panels (f,g) correspond to
$Y = 0.28H$
(just above the canopy tip).
Figure 10(a) presents the scale-energy source term
$\xi$
along with field lines of the flux vector
$\boldsymbol{\varPhi } = (\phi _x, \phi _z, \psi )$
. Repeated indices are omitted here and in the remainder of this section for brevity. The sign of
$\xi$
identifies regions of energy production (
$\xi \gt 0$
) and dissipation (
$\xi \lt 0$
). The flux vector describes the redistribution of energy in scale space (
$r_x$
,
$r_z$
) and the wall-normal direction (
$Y$
); its field lines indicate the average direction of energy transfer, while its divergence,
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varPhi }$
, quantifies the local net flux. A positive divergence implies that fluxes are locally energised, while a negative divergence indicates that fluxes are compensating for local dissipation. Field lines of
$\boldsymbol{\varPhi }$
typically originate in source regions and terminate in sinks, tracing the pathways of energy transfers. However, they do not represent causal trajectories, but rather offer a visual representation of the average energy transfer directions.
The source term
$\xi$
attains its maximum at small scales within the shear layer immediately above the canopy tips, around
$r_x \approx 0$
,
$r_z \in [0.2, 0.25]$
(
$r^{+}_z \in [189, 236]$
), and
$Y \in [0.25, 0.28]$
(
$Y^{+} \in [236, 264]$
). This region coincides with the peak of the production term
$P_{11}$
(not shown), which, owing to flow symmetries, is the sole contributor to scale-energy production. Here, energy is extracted from the mean shear into streamwise velocity fluctuations, while the cross-stream components are sustained through pressure–strain redistribution. This energetic hotspot defines the driving range of scales (Casciola et al. Reference Casciola, Gualtieri, Benzi and Piva2003; Cimarelli et al. Reference Cimarelli, De Angelis and Casciola2013), where strong mean shear (
$\mathrm{d}U/\mathrm{d}y$
) and anisotropy (
$ \left \langle { \delta u \delta v } \right \rangle$
) inject energy into the turbulent field. These production scales coincide with the energy-containing scales identified above the canopy in the previous section, confirming their dynamical relevance. Energy sinks (
$\xi \lt 0$
) are predominantly located within the canopy and at the smallest scales. In contrast to planar channel flows (Gatti et al. Reference Gatti, Chiarini, Cimarelli and Quadrio2020), where the minimum of
$\xi$
occurs at the wall, the strongest dissipation here is found just below the canopy tips. This shift reflects the intense shear and enhanced velocity gradients at the canopy–flow interface, which displace the primary dissipative activity away from the solid boundary. The absence of positive
$\xi$
within the canopy layer indicates that no energy production occurs there; rather, as shown below, the fluctuation structure inside the canopy is largely inherited from motions generated above.
3.2.2. Scale-space energy transfers
We examine the transfer of energy from the production region (
$\xi \gt 0$
) to dissipative regions (
$\xi \lt 0$
) using the scale-space flux vector
$\boldsymbol{\varPhi }$
. The flux organisation reveals a structured redistribution: energy extracted at finite scales within the shear layer above the canopy is conveyed downward into the canopy and upward towards the smallest scales of the outer flow, qualitatively consistent with planar channel turbulence (Cimarelli et al. Reference Cimarelli, De Angelis and Casciola2013).
All flux lines originate from a single topological singularity near
$(r_x, r_z, Y) \approx (0, 0.28, 0.3)$
, i.e.
$(r^{+}_x, r^{+}_z, Y^{+}) \approx (0, 264, 283)$
, coinciding with the previously identified driving range. From this site, two distinct pathways emerge: one directs energy downward into the canopy and ultimately to the wall (see figure 10
b), while the other channels it towards progressively smaller scales in the outer region (see figure 10
c). As dissipative regions are approached, the flux aligns with the local normal direction and becomes increasingly viscous-dominated. In the limit
$r_x, r_z \to 0$
,
$(\phi _x, \phi _z, \psi ) \sim (1,1,0) r$
, supplying small scales to offset wall-parallel dissipation, whereas at the wall
$(\phi _x, \phi _z, \psi ) \sim (0,0,1)Y$
, indicating a wall-normal transfer that compensates dissipation associated with wall-normal gradients.
The downward branch becomes predominantly wall-normal for
$Y \lesssim h$
, revealing a largely scale-independent, top-down transfer (
$\psi \gg \phi _x,\phi _z$
). Energy produced above the canopy is thus injected into the canopy layer, where it is redistributed and ultimately dissipated in the absence of local production. The scales involved (
$r_x \lesssim 0.5$
,
$r_x^+ \lesssim 471$
) are consistent with the KH- and hairpin-type motions identified earlier (§ 3.2.3), indicating that these structures mediate energy transfer between the outer and inner regions. In contrast, the outer branch remains primarily within the overlying flow and follows a reverse-to-direct cascade sequence: energy first shifts towards larger streamwise separations before bending towards the smallest scales, where it feeds dissipation away from the canopy interface. Together, these two pathways delineate a clear energetic partition: a top-down transfer that sustains inner-layer fluctuations, and an outer cascade that redistributes production towards small-scale dissipation aloft.
Energy transfers in scale space can be separated into viscous and turbulent contributions, highlighting the mechanisms governing the cascade (see figure 10
d–g). The viscous flux acts predominantly at the smallest scales, driving a direct cascade along both
$r_x$
and
$r_z$
independently of wall-normal position. Near the origin (
$r_x, r_z \to 0$
), where velocity gradients are strongest, viscous effects dominate and fluxes converge, supplying energy to balance local dissipation. The turbulent flux is dominant across intermediate and larger scales, representing inertia-driven transfer whose magnitude varies with both scale and distance from the wall. Within the canopy (figure 10
$d$
,
$e$
), motions drive an inverse cascade toward larger separations in both directions for
$0 \lt r_x, r_z \lesssim 0.5$
. Just above the canopy (figure 10
$f$
,
$g$
), the process becomes scale-dependent: for
$r_z \lesssim 0.3$
(
$r^{+}_z \lt 283$
), energy first shifts to larger streamwise separations (inverse) before returning to smaller scales (direct), whereas for
$r_z \gtrsim 0.4$
(
$r^{+}_z \gt 377$
) the inverse transfer dominates both directions. Overall, this decomposition shows that viscous fluxes dominate at the smallest scales, while turbulent fluxes control the energy transfer across intermediate and larger scales through combined direct and inverse cascades.
We conclude this section by examining the energetic coupling between the outer flow and the canopy via the average downward-directed flux lines (group I in figure 10). To quantify this exchange locally and across scales, we compute the flux at different streamwise
$X$
and spanwise
$Z$
positions through a circular patch of radius
$r$
in the
$(r_x,r_z)$
plane at the interface
$Y = 0.26 \approx h$
:
\begin{equation} \varPsi (X,Y,Z,r) = \frac {1}{S(r)} \oint _{S(r)} \left ( v^* \delta q^2 + \frac {2}{\rho _f} \delta p \delta v - \frac {\nu }{2} \frac {\partial }{\partial Y} \delta q^2 \right ) \mathrm{d}\varSigma , \quad S(r) = \pi r^2, \end{equation}
where the three terms correspond to turbulent, pressure and viscous contributions. Averaging along the circumference yields similar results. Several scales are considered. The distributions of
$\varPsi$
(see figure 11) indicate, at all scales, a predominant downward transfer of scale energy from the outer flow into the canopy, sustaining motions within it. Upward transfers also occur locally but less frequently. For instance, at
$r = 0.1$
, strong downward events (
$\varPsi \lt -\varPsi _{\textit{rms}}$
) occur with probability 0.111, while upward events (
$\varPsi \gt \varPsi _{\textit{rms}}$
) occur with probability 0.048. The intensity and intermittency of these transfers are largest at smaller scales, reflecting the contribution of intense sweep and ejection events associated with hairpin vortex legs (
$r \approx 0.25$
or
$r^+ \approx 236$
), whereas larger KH-like structures (
$r \approx 0.5$
or
$r^+ \approx 471$
) provide a steadier, more continuous energy supply. Decomposition into turbulent, pressure and viscous components confirms that both downward and upward transfers are dominated by turbulent and pressure-mediated mechanisms, with viscous effects playing a minor role at this wall-normal location. This analysis demonstrates that vertical energy exchange between outer and canopy regions is concentrated in intense, intermittent events rather than uniform transport.
Distribution of the average physical-space flux
$\varPsi$
at
$Y = 0.26H$
(
$Y^{+} = 245$
). The flux is obtained by averaging the instantaneous spatial flux over circular areas of radius
$r$
in the
$r_x$
–
$r_z$
plane, for different values of
$r$
. (a,c) Spatial distribution of the total flux
$\varPsi$
for the scales
$r/H = 0.1, 0.2, 0.3, 0.4, 0.5$
(
$r^{+} \approx 94, 189, 283, 377, 471$
). (b,d) Decomposition of
$\varPsi$
into its constituent components – the turbulent flux
$\varPsi _t$
, pressure flux
$\varPsi _p$
and viscous flux
$\varPsi _v$
.
Figure 11 Long description
Panel A: A probability density function (PDF) graph showing the spatial distribution of the total flux for different values of r/H (0.1, 0.2, 0.3, 0.4, 0.5). The x-axis represents the flux (Ψ) and the y-axis represents the PDF. The graph shows multiple lines representing different r/H values, with the PDF peaking at zero and tapering off symmetrically. Panel B: A PDF graph decomposing the total flux into its constituent components: turbulent flux (Ψ_t), pressure flux (Ψ_p), and viscous flux (Ψ_v). The x-axis represents the flux (Ψ) and the y-axis represents the PDF. The graph shows three lines representing the different flux components, each peaking at zero and tapering off symmetrically. Panel C: A PDF graph showing the spatial distribution of the total flux normalized by its root mean square value (Ψ/Ψ_rms) for different values of r/H (0.1, 0.2, 0.3, 0.4, 0.5). The x-axis represents the normalized flux (Ψ/Ψ_rms) and the y-axis represents the PDF. The graph shows multiple lines representing different r/H values, with the PDF peaking at zero and tapering off symmetrically. Panel D: A PDF graph decomposing the normalized total flux into its constituent components: turbulent flux (Ψ_t), pressure flux (Ψ_p), and viscous flux (Ψ_v). The x-axis represents the normalized flux (Ψ/Ψ_rms) and the y-axis represents the PDF. The graph shows three lines representing the different flux components, each peaking at zero and tapering off symmetrically.
3.2.3. The role of pressure–strain in shaping the canopy-layer fluctuation structure
(a–c) Pressure–strain term
$\varPi _{ij}$
in the
$r_y = 0$
plane, zoomed inside the canopy layer for
$Y \lt 0.2$
. White lines indicate contours of
$\varPi _{ij}$
at various levels, with the black line marking the
$\varPi _{ij} = 0$
contour, and are used to identify the local minima/maxima. (d) The single-point pressure–strain terms as a function of the wall-normal coordinate
$Y$
. (e) An example of energy redistribution within the canopy layer driven by pressure–strain. The isosurface corresponds to negative Fourier-filtered pressure fluctuations at
$\tilde {p} = -0.35$
, superimposed with vectors of the Fourier-filtered velocity field
$\boldsymbol{\tilde {u}}$
, coloured by the intensity of the vertical component
$\tilde {v}$
.
We now examine the mechanism shaping the structural organisation of velocity fluctuations within the canopy layer (§ 3.1). In this region,
$\xi \lt 0$
, indicating that the local flow is inactive in producing turbulent energy. Instead, the scale-dependent structure of the fluctuations is inherited from motions above the canopy, primarily through sweep events, and redistributed locally by the pressure–strain term,
$\varPi _{ij}$
.
Figure 12 shows the scale-space distribution of
$\varPi _{ij}$
within the canopy (
$Y \leqslant 0.2$
). Across all scales and positions,
$\varPi _{22} \lt 0$
while
$\varPi _{11} \gt 0$
and
$\varPi _{33} \gt 0$
. The no-penetration boundary condition enforces an effectively two-dimensional near-wall flow, so vertical velocity fluctuations are reoriented towards wall-parallel directions. Compared with canonical planar channel flows (Gatti et al. Reference Gatti, Chiarini, Cimarelli and Quadrio2020), the region of
$\varPi _{22} \lt 0$
extends further from the wall, reflecting enhanced suppression of vertical motions by the canopy filaments. Within the canopy,
$\varPi _{22}$
exhibits two localised minima at
$(Y, r_x, r_z) \approx (0.013, 0, 0.65{-}0.8)$
and
$(0.013, 0.5{-}0.65, 0)$
, coinciding with positive peaks in
$\varPi _{11}$
and
$\varPi _{33}$
at
$(Y, r_x, r_z) \approx (0.012,0.45{-}0.55,0)$
and
$(0.013{-}0.025,0,0.4{-}0.65)$
. These locations correspond to the dominant flow scales identified near the canopy tip (see § 3.1 and figure 6).
Despite attenuation by the filaments, the canopy-layer fluctuations retain a scale-dependent structure dictated by the outer flow and mediated by pressure–strain redistribution. Streamwise fluctuations (
$u$
) reflect the characteristic length scales of KH rolls, while spanwise fluctuations (
$w$
) preserve the spanwise scales associated with the legs of hairpin-like vortices. This provides a direct mechanistic link between outer-layer dynamics and the observed organisation within the canopy, confirming that pressure–strain-mediated redistribution governs the scale-specific imprint of the outer flow. An instantaneous snapshot (figure 12
e) illustrates the mechanism: large-scale sweeps generated by KH rolls inject high-momentum fluid into the canopy. Interaction with the filaments and wall redirects these vertical motions via pressure–strain, enhancing wall-parallel
$u$
fluctuations at scales consistent with the KH rolls. A similar pattern occurs for
$w$
fluctuations at
$r_z \ne 0$
, highlighting the role of hairpin legs in mediating spanwise momentum redistribution.
3.3. Flexibility effects
Influence of filament flexibility on outer-layer flow scales. Results are shown for three Cauchy numbers:
$\textit{Ca} = 0$
,
$\textit{Ca} = 25$
and
$\textit{Ca} = 50$
. (a) Correlation coefficient
$\rho _{22}$
as a function of
$r_x$
in the
$r_y = r_z = 0$
plane, evaluated at the wall-normal location where
$\langle \delta v \delta v \rangle$
peaks. (b) Same as (a), but in the
$r_x = r_y = 0$
plane. (c) Same as (b), but for
$\rho _{11}$
and
$\langle \delta u \delta u \rangle$
. (d) Turbulent production term
$P_{11}$
as a function of
$r_z$
in the
$r_x = r_y = 0$
plane, evaluated at the wall-normal location corresponding to the peak of
$\langle \delta u \delta u \rangle$
. In all panels, insets show the corresponding quantities in viscous units.
Figure 13 Long description
Panel A: A line graph shows the correlation coefficient as a function of r_x in the x-z plane, evaluated at the wall-normal location where ρ_22 peaks. The graph includes three lines representing different Cauchy numbers: Ca = 0 (blue), Ca = 25 (magenta), and Ca = 50 (yellow). The inset shows the corresponding quantities in viscous units. Panel B: A line graph shows the correlation coefficient as a function of r_z in the x-z plane. The graph includes the same three lines representing different Cauchy numbers. The inset shows the corresponding quantities in viscous units. Panel C: A line graph shows the correlation coefficient ρ_11 as a function of r_z in the x-z plane. The graph includes the same three lines representing different Cauchy numbers. The inset shows the corresponding quantities in viscous units. Panel D: A line graph shows the turbulent production term P_11 as a function of r_z in the x-z plane, evaluated at the wall-normal location corresponding to the peak of ρ_11. The graph includes the same three lines representing different Cauchy numbers. The inset shows the corresponding quantities in viscous units.
We now assess the role of filament flexibility across
$\textit{Ca} \in [0,50]$
. While the dominant physical mechanisms and key flow structures remain broadly consistent across canopy configurations, increasing
$\textit{Ca}$
introduces systematic quantitative differences. As shown in figure 2(b), higher
$\textit{Ca}$
enhances turbulent kinetic energy within the canopy due to increased filament compliance, while energy above the canopy is slightly reduced, reflecting a weaker and more diffuse shear layer.
To quantify these effects, we focus on two aspects. First, we assess the coherence of the flow structures as
$\textit{Ca}$
varies, revealing that more flexible canopies tend to partially decorrelate the KH-like rolls and hairpin structures, particularly near the canopy tips. Second, we examine how the canopy modulates vertical energy transfer from the outer flow into the interior. Increased filament flexibility effectively creates a partially shielding layer, attenuating the penetration of outer-layer motions and modifying the scale-dependent energy fluxes within the canopy. This analysis highlights that while the flow retains its qualitative organisation, the quantitative redistribution of turbulent energy across scales is sensitive to canopy compliance, with implications for both momentum transport and the imprint of coherent structures inside the canopy.
The influence of
$\textit{Ca}$
on outer-layer flow structures is quantified in figure 13, which reports the two-point correlation coefficients
$\rho _{11}$
and
$\rho _{22}$
as functions of
$r_x$
and
$r_z$
, evaluated at the wall-normal locations corresponding to the peaks of the respective structure functions. With increasing
$\textit{Ca}$
, the peak of
$\rho _{22}(r_x)$
(figure 13
a), associated with the KH rolls, decreases in magnitude and shifts towards larger
$r_x$
and lower wall-normal positions. This reflects a weakening and broadening of the shear layer, producing KH rolls that are less coherent and more diffuse. For
$\textit{Ca} = 50$
,
$\rho _{22}(r_x)$
becomes only marginally negative, indicating near-complete suppression of KH activity. Figure 13(b,c) shows
$\rho _{22}$
and
$\rho _{11}$
as functions of
$r_z$
, capturing the statistical footprints of streamwise-oriented vortices and the LSS and HSS. The coherence and intensity of the LSS/HSS decrease with increasing
$\textit{Ca}$
(consistent with Sundin & Bagheri (Reference Sundin and Bagheri2019)), with a slight reduction of the spanwise scale in both outer and viscous units. In contrast, the signature of streamwise-oriented vortices is only marginally affected by canopy flexibility, with the peak of
$\rho _{22}$
slightly moving towards smaller spanwise scales as
$\textit{Ca}$
increases. The attenuation of the LSS/HSS is linked to a reduction in turbulent production
$P_{11}$
, resulting from the weakened shear layer (i.e. smaller
$\mathrm{d}U/\mathrm{d}y$
; see figure 2
a). In the most flexible case (
$\textit{Ca} = 50$
), the secondary peak in
$\langle \delta u \delta u \rangle$
in the
$r_x = r_y = 0$
plane, prominent in the rigid canopy (figure 5), disappears, further confirming the near-complete suppression of KH rolls.
Influence of filament flexibility on canopy-layer flow scales. Results are shown for three Cauchy numbers:
$\textit{Ca} = 0$
,
$\textit{Ca} = 25$
and
$\textit{Ca} = 50$
. (a) Correlation coefficient
$\rho _{11}$
as a function of
$r_x$
in the
$r_y = r_z = 0$
plane, evaluated at
$Y = 0.03H$
. (b) Same as (a), but for
$\rho _{33}$
in the
$r_x = r_y = 0$
plane. (c) Pressure–strain term
$\varPi _{11}$
as a function of
$r_x$
in the
$r_z = r_y = 0$
plane, evaluated at
$Y = 0.03H$
. (d) Same as (c), but for
$\varPi _{33}$
in the
$r_y = r_x = 0$
plane.
Within the canopy layer, the flow remains dominated by streamwise and spanwise fluctuations for all
$\textit{Ca}$
considered, whose organisation reflects the imprint of the overlying motions near the canopy tips (see figure 12). Figure 14(a,b) shows
$\rho _{11}$
and
$\rho _{33}$
as functions of
$r_x$
and
$r_z$
at
$Y = 0.03$
, with similar trends for small variations in
$Y$
. Fluctuation intensity increases with
$\textit{Ca}$
across all scales, accompanied by stronger pressure–strain redistribution (figure 2;
$\varPi _{11}$
and
$\varPi _{33}$
in figure 14
c,d), indicating enhanced reorientation of wall-normal motions into wall-parallel components, likely due to stronger vertical confinement. At the same time, the spatial coherence of the fluctuations decreases with flexibility, reflecting a loss of structural organisation consistent with the attenuation of outer-layer motions.
Influence of filament flexibility on the distribution of the average physical-space flux
$\varPsi$
, evaluated at the mean canopy tip location
$Y = \langle y_{{tip}} \rangle _f$
for each case. The flux is computed as described in § 3.2. Results are shown for three Cauchy numbers:
$\textit{Ca} = 0$
,
$\textit{Ca} = 25$
and
$\textit{Ca} = 50$
. (a,c) Spatial distribution of the total flux
$\varPsi$
at scale
$r/H = 0.1$
. (b,d) Spatial distribution of the total flux
$\varPsi$
at scale
$r/H = 0.5$
.
Figure 15 Long description
Panel A: A line graph shows the spatial distribution of the total flux at scale for three Cauchy numbers: 0, 25, and 50. The x-axis represents the flux (Ψ) ranging from -1.0 to 1.0, and the y-axis represents the probability density function (PDF) on a logarithmic scale from 10^-5 to 10^0. Panel B: Another line graph displays the spatial distribution of the total flux at scale for the same Cauchy numbers. The x-axis represents the flux (Ψ) ranging from -1.0 to 1.0, and the y-axis represents the PDF on a logarithmic scale from 10^-5 to 10^0. Panel C: A line graph illustrates the spatial distribution of the total flux at scale normalized by the root mean square (Ψ/Ψrms) for the three Cauchy numbers. The x-axis ranges from -50 to 50, and the y-axis represents the PDF on a logarithmic scale from 10^-5 to 10^0. Panel D: The final line graph shows the spatial distribution of the total flux at scale normalized by the root mean square (Ψ/Ψrms) for the same Cauchy numbers. The x-axis ranges from -50 to 50, and the y-axis represents the PDF on a logarithmic scale from 10^-5 to 10^0.
We conclude this section by examining how filament flexibility affects the energetic coupling between the outer flow and the canopy. Figure 15 reports the scale-averaged flux
$\varPsi$
at two representative scales,
$r = 0.1$
and
$r = 0.5$
, evaluated at the mean filament-tip height. For all
$\textit{Ca}$
, the distributions remain left-skewed, indicating a net top-to-bottom energy transfer across scales (see § 3.2.1). Increasing
$\textit{Ca}$
reduces the magnitude of the distribution tails (figure 15
a,b), reflecting a weaker shear layer and a loss of coherence in sweep–ejection structures. Energy transfer is strongest at smaller scales, where the influence of the flexibility is most pronounced. When normalised by the root-mean-square value (figure 15
c,d), the distributions collapse across all
$\textit{Ca}$
, indicating that the relative intermittency of energy transfer remains largely unaffected by filament flexibility.
4. Conclusions
Schematic representation of the mechanism sustaining velocity fluctuations in canopy turbulence.
This study has addressed the central question of how turbulence is organised and sustained in dense, submerged canopy flows, and how energy is exchanged across scales and between the canopy layer and the overlying shear flow. Using the DNS database of Monti et al. (Reference Monti, Olivieri and Rosti2023) and the AGKE (Gatti et al. Reference Gatti, Chiarini, Cimarelli and Quadrio2020), we have provided a space- and scale-resolved characterisation of production, redistribution and inter-scale transfer for both rigid and flexible canopies at
$\textit{Re}_b = 5000$
and Cauchy numbers
$\textit{Ca} = 0$
,
$25$
and
$50$
.
Figure 16 summarises the main energy pathways across scales and flow regions. Energy production is localised at the interfacial shear layer, over a narrow range of streamwise and spanwise scales aligned with the dominant coherent structures in the outer layer, such as KH-like rolls, hairpin vortices and streaks. Within the canopy, fluctuations are sustained non-locally, primarily through inter-scale transfer and pressure–strain redistribution, inheriting their organisation and characteristic scales from the outer-layer motions. Sweep events induced by the outer-layer structures play a key role in transferring energy downward, while pressure–strain interactions redistribute it among velocity components, shaping the anisotropy of canopy-layer fluctuations.
The energy exchange across the canopy interface is asymmetric but bidirectional. While the net transfer is from the outer layer into the canopy, intermittent reverse fluxes occur at all scales, with the strongest events dominated by small-scale motions. This interplay establishes a dynamic coupling between inner and outer regions, linking coherent outer-layer structures to the organisation of canopy turbulence. Vertical energy exchange at the canopy interface is therefore concentrated in intense, intermittent events rather than uniform transport. In contrast, velocity fluctuations in the bulk flow are sustained by upward transfers originating from the shear layer above the canopy tip, combining both inverse and direct cascades towards the smallest scales.
Canopy flexibility modifies both outer- and inner-layer dynamics. Increasing compliance weakens the coherence of outer-layer structures, diffuses the shear layer and reduces the efficiency of inter-layer energy exchanges. Within the canopy, although the overall energy increases, the coherence of structures inherited from the motions above is diminished, and the scale-dependent coupling with the outer flow is attenuated, leading to weaker energy exchanges across the canopy interface.
Overall, this study provides a unified energetic interpretation that connects coherent structures to production, redistribution and inter-scale transfer mechanisms. Future work could extend this analysis to broader parameter ranges, for example by exploring sparse canopies, where large-scale structures penetrate deeper into the canopy layer (Sharma & García-Mayoral Reference Sharma and García-Mayoral2020b ), varying canopy height or artificially constraining filaments to be rigid to isolate the effect of reconfiguration (Foggi Rota et al. Reference Foggi Rota, Monti, Olivieri and Rosti2024). It would also be of interest to investigate non-uniform or patchy canopy arrangements (Löhrer & Fröhlich Reference Löhrer and Fröhlich2023; Foggi Rota et al. Reference Foggi Rota, Tressoldi, Avallone and Rosti2025), where added inhomogeneity is expected to significantly alter the scale-resolved energy transfer mechanisms identified in this study.
Acknowledgements
The authors acknowledge the computer time provided by the Scientific Computing and Data Analysis section of the Core Facilities at OIST and the computational resources offered by the HPCI System Research Project under grant numbers hp210025, hp220402, hp240006 and hp250035.
Funding
The research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding to M.E.R. from the Cabinet Office, Government of Japan. M.E.R. acknowledges funding from the Japan Society for the Promotion of Science (JSPS, grant numbers 24K00810 and 24K17210). R.B. acknowledges partial funding from Politecnico di Milano.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Fluid velocity on the filaments
In this appendix, we show that the immersed-boundary method employed in the present study (Yu Reference Yu2005; Huang et al. Reference Huang, Shin and Sung2007) rigorously enforces the no-slip and no-penetration boundary conditions at the Lagrangian points representing the canopy filaments. To provide the most stringent test, we consider the configuration featuring rigid filaments, as this induces the strongest hydrodynamic back-reaction on the fluid. We extract the velocity field from an instantaneous snapshot along a spanwise line (parallel to the
$z$
axis) passing directly through a representative stem located at
$(x/H, z/H) = (0.11, 0.36)$
. As shown in figure 17, all three velocity components (
$u, v, w$
) vanish exactly at the coordinate of the Lagrangian point, confirming that the required kinematic constraints are strictly satisfied. The Lagrangian force responsible for enforcing this condition is projected onto the underlying Eulerian fluid grid over a three-point compact support via convolution with a regularised delta function (Roma et al. Reference Roma, Peskin and Berger1999).
Instantaneous velocity components extracted along a spanwise line passing through the representative rigid stem located at
$(x/H, z/H) = (0.11, 0.36)$
. All three velocity components vanish exactly at the filament coordinate. Large markers indicate the Eulerian grid points within the compact support of the immersed-boundary-method forcing.
Reb
Lx
Lz
Nx
Ny
Nz
nx
nz
λ
⟨dp/dx⟩
⟨ytip⟩f
⟨⋅⟩f
U
k
⟨uu⟩
⟨vv⟩
⟨ww⟩
Ca=0
Ub
H
⟨ytip⟩f/H=0.25,0.19,0.135
Ca=0,25,50
⟨⋅⟩f
Ca=0
u=±0.6
p~=−0.06
u
⟨δuiδuj⟩
ry=0
Ca=0
⟨δuδu⟩
⟨δvδv⟩
⟨δwδw⟩
⟨δuiδuj⟩
ry=rz=0
rx=ry=0
Rij
⟨δuiδuj⟩
Rij
⟨δuδu⟩
R11
(Y,rz)
Ca=0
R11=0
Ca=0
⟨δuiδuj⟩
ry=0
Y<0.2
⟨δuiδuj⟩
⟨δuiδuj⟩
Y
Ca=0
Ca=0
p~=−0.04
v~
u~v~<0
y≈h
Q
Q~=1.5
ω~z
p^=−0.065
v^
Ca=0
Q^=1
ω^x
u^2D=(v^,w^)
v^
ξ
ry=0
Ca=0
Φ=(ϕx,ϕz,ψ)
∇⋅Φ
ξ=0.9ξmax
rx
rz
Y
Φ
s
s=1/smax∫0smaxds
ds=drx2+drz2+dY2
∇⋅Φ<0
∇⋅Φ>0
Φt
Φv
ry=0
Y=0.15H
Y=0.28H
Ψ
Y=0.26H
Y+=245
r
rx
rz
r
Ψ
r/H=0.1,0.2,0.3,0.4,0.5
r+≈94,189,283,377,471
Ψ
Ψt
Ψp
Ψv
Πij
ry=0
Y<0.2
Πij
Πij=0
Y
p~=−0.35
u~
v~
Ca=0
Ca=25
Ca=50
ρ22
rx
ry=rz=0
⟨δvδv⟩
rx=ry=0
ρ11
⟨δuδu⟩
P11
rz
rx=ry=0
⟨δuδu⟩
Ca=0
Ca=25
Ca=50
ρ11
rx
ry=rz=0
Y=0.03H
ρ33
rx=ry=0
Π11
rx
rz=ry=0
Y=0.03H
Π33
ry=rx=0
Ψ
Y=⟨ytip⟩f
Ca=0
Ca=25
Ca=50
Ψ
r/H=0.1
Ψ
r/H=0.5
(x/H,z/H)=(0.11,0.36)
