3.1. Mean flow characteristics

We begin by examining the effect of in-plane solidity on the mean velocity profile. Figure 3(a) presents the dimensionless mean velocity distributions normalised by the open channel height ($H$) and the bulk velocity ($U_b=({1}/{H})\int _0^H\langle U \rangle \,{{\rm d} y}$, where $\langle {\cdot } \rangle$ represents the triple averaging operator encompassing time and homogeneous directions).

Figure 3. (a) Mean velocity profiles made non-dimensional with the canopy height $h$ and the bulk velocity $U_b$. See table 1 for symbols. The red and blue symbols represent the internal inflection point and the virtual origin, respectively. (b) Comparison with the mean profiles inside the canopy obtained by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) for the same solidity values: blue solid line is the present HD case ($\lambda =0.56$); red solid line is the present DE case ($\lambda =0.32$); blue dashed line is Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) $\lambda =0.56$ case and red dashed line is Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) $\lambda =0.32$ case.

Distinct variations in the mean velocity distribution are evident for different solidity values $\lambda$, particularly within the canopy region ($\kern0.7pt y/H \in [0,0.1]$). As the solidity increases, the mean velocity profiles exhibit increased retardation with a progressively enhanced concave shape. The sparse (SP) case demonstrates a minor slowdown, while the TR, MD, DE cases (i.e. cases with increasing solidity) exhibit mild retardation. In contrast, the densest (HD) configuration demonstrates a substantial decrease in mean velocity accompanied by significant curvature enhancement.

Comparing the present HD profile with the one obtained by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) at the same $\lambda$ value (see figure 3b), clear differences in shape are evident. Specifically, the velocity gradient in the wall region is much milder in the present case, with a smoother change of curvature when moving away from the wall. The strong gradient in the wall region and the appearance of a kink in the velocity profile, just below the inner inflection point, reported by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020), is absent in the present mean velocity profile.

The presence of the kink was attributed to the strong inrushes (jet-like) of wall-normal velocity (Banyassady & Piomelli Reference Banyassady and Piomelli2015), deflected in the $x$ and $z$ directions due to continuity and wall impermeability. Although the two dense cases share the same solidity $\lambda$, the absence or presence of the kink suggests potential differences in the flow dynamics within the canopy region.

Furthermore, it is interesting to note that the mean velocity profile of the present HD case is similar to those reported in the context of turbulent wall flows bounded by porous surfaces (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006; Rosti, Cortelezzi & Quadrio Reference Rosti, Cortelezzi and Quadrio2015; Kuwata & Suga Reference Kuwata and Suga2017).

When considering a non sparse canopy flow regime (i.e. $\lambda >0.04$), the mean velocity profile features three notable points within the filamentous layer (i.e. for $y\le h$ as shown in figure 3a): two inflection points and the location of the virtual origin of the outer flow (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004; Nepf Reference Nepf2012; Monti et al. Reference Monti, Omidyeganeh and Pinelli2019). The latter can be interpreted as the location of the virtual wall seen by the outer turbulent boundary layer above the canopy. The concept of a shifted boundary layer is commonly proposed in the framework of boundary layers developing over textured surfaces (Jiménez Reference Jiménez2004; Bottaro Reference Bottaro2019). In our case, to determine the location of the virtual origin ($\kern0.7pt y_{vo}$), it is assumed that the outer mean velocity profile follows the canonical logarithmic law (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020):

(3.1)\begin{equation} \frac{\langle U \rangle}{u_{\tau \,out}}=\frac{1}{\kappa}\log\left(\frac{y_{out}\,u_{\tau \,out}}{\nu} \right) +B, \end{equation}

where $y_{out}=y-y_{vo}$ is the shifted wall-normal coordinate (i.e. $y$ measured from the virtual origin), $\kappa$ is the von Kármán constant and $B$ is a shift accounting for the nature of the wall roughness (Jiménez Reference Jiménez2004). In the present work, the friction velocity $u_{\tau,out}$ in the above equation is computed from the total stress at the location of the virtual origin as the sum of the viscous stresses and the Reynolds shear stress:

(3.2)\begin{equation} u_{\tau,out}=\sqrt{\frac{\tau(\kern0.7pt y_{vo})}{\rho}}=\sqrt{\nu\frac{\partial \langle U \rangle }{\partial y}_{vo}-\langle u'v'\rangle_{vo}}. \end{equation}

By ensuring the alignment of the von Kármán constant $\kappa$ (here assumed to be $\kappa =0.41$), the wall-normal position of the virtual origin for each canopy configuration can be determined if the mean velocity and total stress profile are known. The intercept value $B$ in the logarithmic law is not an unknown, instead, it is a result obtained from solving (3.1) for $y_{vo}$. Additional information regarding this process, as well as the sensitivity of the virtual origin's location to the chosen value of $\kappa$, can be found in the study by Monti et al. (Reference Monti, Omidyeganeh and Pinelli2019).

Figure 3(a) reveals the presence of two inflection points on each mean velocity profile: an outer inflection point located at the top of the canopy and an internal inflection point near the wall. The outer inflection point is generated by the abrupt discontinuity of the drag force exerted by the canopy at its tip. The internal inflection point results from the merging of the velocity profile close to the wall with the mean velocity distribution in the region closer to the canopy tip (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020).

The position of the inner inflection point can be determined by finding the location where the second derivative of the mean velocity is zero. This can be accomplished by considering the mean momentum equation in the streamwise direction:

(3.3)\begin{equation} \frac{\partial P}{\partial x}=\frac{1}{Re_b}\frac{\partial^2 \langle U \rangle }{\partial y^2}-\frac{\partial \langle u'v'\rangle}{\partial y}-\langle F_D\rangle. \end{equation}

In (3.3), the first, second and third terms represent the mean streamwise pressure gradient, the viscous force and the Reynolds shear stress, respectively. The final term accounts for the overall mean drag imposed by the canopy filaments and, since the filaments are mounted normally to the bottom wall, its contribution is mainly due to the form drag of the cylinders, with a weak contribution from the viscous shear stress.

Figure 4 illustrates the mean locations of the notable points of the mean velocity profiles, namely the inner inflection point (red horizontal line), virtual origin (black horizontal line) and outer inflection point (blue horizontal line), as a function of the solidity $\lambda$ for the five canopy configurations considered in this work. Note that when the virtual origin moves below the inner inflection point, no overlapping region can be defined any longer.

Figure 4. Mean locations of the notable points of the mean velocity profiles: inner inflection point (red horizontal line); virtual origin (black horizontal line) and the outer inflection point (blue horizontal line) as a function of the solidity $\lambda$. Note that the canopy becomes sparser (smaller $\lambda$) when increasing the mean filament spacing $\Delta S$.

In the case of sparse canopies, characterised by $\Delta S /h \gg 1$ (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004; Sharma & García-Mayoral Reference Sharma and García-Mayoral2020a), as $\lambda$ is decreased, the virtual origin moves towards the canopy bed (i.e. $y_{vo}/H\rightarrow 0$) and the internal inflection point collapses onto the canopy tip. Conversely, as $\lambda$ is increased, the location of the virtual origin starts moving away from the canopy bed and the inner inflection point moves towards it. The signed distance between these two points has a zero in the range $0.14<\lambda <0.25$. The value of $\lambda$ at which the two points collapse into a single location can be used as a criterion for predicting the onset of the dense canopy flow regime (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020).

In dense configurations, the virtual origin is located in the upper half of the canopy (i.e. $y_{vo}/h>0.5$), while the inner inflection point is in the bottom half of the canopy (i.e. $y_{vo}/h<0.5$). An increase in the width of the overlap region suggests a decoupling between the inner and outer layers (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020; Sharma & García-Mayoral Reference Sharma and García-Mayoral2020b). In the densest, HD case, the virtual origin and internal inflection points are located near the canopy tip and bed, respectively. In this situation, the outer flow behaves like a shear flow over a rough wall, while the innermost flow develops a local boundary layer near the bed.

Figure 5(a,b) show the locations of the inner inflection point and virtual origin (the latter measured from the canopy tip) as a function of $\lambda$. These figures, which incorporate data from Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020), demonstrate that the parameter $h-y_{vo}$ and the location of the internal inflection point saturate to a constant value as $\lambda$ increases (denser canopies). This behaviour can be better understood by considering a scenario where the filament spacing is kept constant (i.e. $\Delta S/H={\rm \pi} /24$) while systematically increasing the canopy height. In this case, the inner and outer layers gradually decouple until reaching conditions where the canopy height (or the solidity $\lambda$ in this case) becomes irrelevant to the flow behaviour.

Figure 5. (a) Mean locations of the inner inflection point as a function of the solidity $\lambda$. (b) Mean locations of the virtual origin measured from the canopy tip as a function of $\lambda$. Both parameters are made non-dimensional using the half-channel height $H$. The solid lines are best fits for the current data, while the dashed lines fit the results reported by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). Symbols meanings are in table 1.

However, in the current framework where $\lambda$ is increased by considering progressively denser canopies of the same height, no clear saturation is observed. The latter will eventually occur when extremely packed filaments are considered, i.e. almost solid wall as $\lambda \rightarrow \infty$.

Further information on the locations of the notable points can be extracted by looking at their scaling with the canopy geometric parameters. Figure 6 shows the locations of the internal inflection points $y_{{in}}$ of the mean velocity profiles, scaled with $\Delta S$ and $h$, plotted versus $\lambda$. Figure 6(a) demonstrates a collapse between the cases sharing the same $\lambda$ value, regardless of whether it is obtained by changing $\Delta S$ or $h$. As $\Delta S$ approaches zero (i.e. increasing $\lambda$), the wall becomes solid and $y_{{in}}/h$ approaches zero as well, resulting in a linear decreasing trend. However, when $h$ becomes very large, the behaviour of the inner region flow becomes independent of $h$ (or $\lambda$) and $y/h$ reaches an asymptotic value.

Figure 6. (a) Mean locations of the inner inflection point made dimensionless using $\Delta S$ as a function of the solidity $\lambda$. (b) Mean locations of the virtual origin made dimensionless with the canopy height $h$ as a function of $\lambda$. The solid lines are best fits for the current data, while the dashed line fits the results reported by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). Symbols meanings are in table 1.

Figure 7 provides insights into the mean location of the virtual origin measured from the canopy tip as a function of $\lambda$. Unlike the thickness of the inner layer, which scales with $\Delta S$, the thickness of the outer layer scales with the canopy height $h$. The open symbols in figure 7(a) indicate that the outer thickness remains roughly the same as the filament average spacing, indicating that it is $\Delta S$ that determines the size of the eddies able to penetrate the canopy (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020).

Figure 7. Mean location of the virtual origin measured from the canopy tip, as a function of solidity ($\lambda$). (a) Scaling with $\Delta S$ and (b) scaling with the canopy height $h$. The solid line represents a polynomial fit passing through the virtual origin positions of the full canopy database (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020).

However, when $\Delta S$ is kept constant, the penetration is governed by $h$, and asymptotically saturates to unity as the canopy height increases towards the densest regime. Figure 7(b) shows that when $(h-y_{{vo}})$ is scaled with $h$, a collapse is observed irrespective of how the value of $\lambda$ is obtained. This result is consistent with the interpretation of $(h-y_{{vo}})/h$ as the fraction of the canopy that the outer eddies can penetrate. In summary, $\Delta S$ determines the size of the eddies that can enter the canopy, while $h$ determines the depth of penetration.

3.2. Notable points locations: heuristic model

We propose a simple geometric model to explain the variations in the position of the virtual origin of the mean velocity profile (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). Additionally, we discuss the influence of the mean filament spacing on the upwash/downwash events produced inside the canopy by outer turbulent eddies and their role in determining the transition between different canopy flow regimes (Nepf Reference Nepf2012). By keeping the canopy height $h$ constant, the variations in $\Delta S$ can be interpreted as modulations of the size of the high-pass filter acting on the outer flow structures. In particular, as illustrated in figure 8, it is the ratio $(\Delta S - d)/h$ that determines the size of the eddies capable of penetrating the canopy and their depth of penetration. For slender canopies ($d/h \ll 1$), such as the cases under consideration in this work, the ratio $\Delta S /h$ is as a sufficient indicator.

Figure 8. Sketch illustrating the largest possible vortex size that can penetrate the canopy from the outer layer and the corresponding location of the virtual origin. An individual vortex is depicted as an idealised eddy with diameter $\Delta S -d$.

In the case of $\Delta S/h > 1$, when the width of the canopy voids exceeds the characteristic eddy size (as shown in figure 8a), the mean filament spacing does not impose any constraints on the eddy size that can penetrate the canopy. In this situation, the mean location of the virtual origin, $y_{{vo}}$, is determined by the average distance between the eddy core and the impermeable bottom wall (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020).

Conversely, in a dense canopy regime, when $\Delta S/h < 1$, only eddies with sizes smaller than $\Delta S$ can fit within the void regions of the canopy (see figure 8c). If an external eddy exceeds $\Delta S$, its interaction with the filaments leads to a local deformation, resulting in a penetration length smaller than $\Delta S$. In this situation, the virtual origin position, seen by the outer flow, moves closer to the canopy tip at a distance of approximately $\Delta S$ from it. The mean location of $y_{{vo}}$ would scale directly with $\Delta S$ (see figure 9), which supports the hypothesis proposed by MacDonald et al. (Reference MacDonald, Ooi, García-Mayoral, Hutchins and Chung2018) for turbulent flows over spanwise aligned bars. The transitional state between the dense and sparse asymptotic conditions occurs when the canopy aspect ratio approaches unity ($\Delta S/h \approx 1$). In this regime, the void regions can be approximated as square cavities accommodating eddies with a size of $\approx \Delta S$, resembling turbulent flows over d-type rough surfaces (Perry, Schofield & Joubert Reference Perry, Schofield and Joubert1969; Leonardi, Orlandi & Antonia Reference Leonardi, Orlandi and Antonia2007; MacDonald et al. Reference MacDonald, Chan, Chung, Hutchins and Ooi2016). The mean location of the virtual origin in this case is positioned at the midpoint of the canopy height (as visible from figure 9), supporting the proposed transition criterion.

Figure 9. Mean location of the virtual origin scaled with the mean filament spacing $\Delta S$, as a function of the aspect ratio $\Delta S/h$. See table 1 for symbols meaning.

Figure 10 provides a qualitative assessment of the transition criterion, showing the contours of streamwise and wall-normal vorticity fluctuations on a $y$–$z$ cross-section of the channel for different $\lambda$ values. In the sparser regime, the outer flow fully penetrates the canopy layer, while in the denser regime, the penetration depth of the outer vortical structures is influenced by $\Delta S$. The localised fragmentation of large-scale vortical structures at the canopy interface in the densest cases aligns with the model predictions. By establishing the criterion for the onset of the dense regime ($\Delta S - d \approx h$), we can provide a precise estimate of the mean filament spacing and the corresponding solidity value ($\lambda$) that leads to the canopy flow transition. Considering the specific geometric parameters of the canopy used in this study ($h/H = 0.1$ and $d/h = 0.24$), we obtain $\Delta S/H = 0.124$, corresponding to $\lambda =(hd)/\Delta S^2 = (H/\Delta S)^2 (d/h)( h/H) \approx 0.15$ or $\Delta S/h =1+d/h\approx 1.24$. The transitional $\lambda$ value of $0.15$ closely aligns with previous findings by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). It is also noticed that a similar threshold value of $\lambda \approx 0.15$ is commonly accepted for distinguishing sparse from dense regimes in k-type rough surfaces (Schlichting Reference Schlichting1937; MacDonald et al. Reference MacDonald, Chan, Chung, Hutchins and Ooi2016).

Figure 10. Contours of the instantaneous vorticity fluctuations extracted from a $y$–$z$ cross-plane: (a,c,e,g,i) $\omega _x'$ (i.e. streamwise component); (b,d,f,h,j) $\omega _y'$ (i.e. wall-normal component). Each row, from top to bottom, represents an increasing solidity value (see table 1). The colour code, from red to blue, represents positive to negative values. The vorticity fluctuations have been made non-dimensional with the outer flow variables (i.e. $H$ and $U_b$), and the corresponding dimensionless vorticity fluctuations $\omega '_x H/U_b$ and $\omega '_y H/U_b$ are both within the range $[-8,8]$.

3.3. Statistical characterisation of the flow

We examine the mean velocity profiles by separately considering the region within the filamentous layer from the bed to the virtual origin and the outer region above $y_{{vo}}$. Figure 11(a) presents the mean velocity distribution in a semi-logarithmic scale for all the solidity values. The profiles are normalised using two friction velocities: $u_{\tau,{in}}=\sqrt {\tau _w/\rho }$ at the impermeable wall ($\kern0.7pt y=0$) and $u_{\tau,{out}}=\sqrt {\tau (\kern0.7pt y_{{vo}})/\rho }$ at the virtual origin ($\kern0.7pt y=y_{{vo}}$), where $\tau _w$ is the wall shear stress at the bed and $\tau (\kern0.7pt y_{{vo}})$ is the total stress at the virtual origin. The figure shows that regardless of the solidity value, the mean velocity profiles collapse within the first few wall units $\tilde {y}^+$ (see figure caption for its definition), indicating the dominance of wall friction over the drag from filaments, even for denser canopies. As we move away from the wall, the modulation imposed by the filaments and the effects of the external flow start to play a role in the mean velocity distribution.

Figure 11. (a) Mean velocity profiles normalised with the two friction velocities $u_{\tau,in}=\sqrt {\tau _w/\rho }$ and $u_{\tau,out}= \sqrt {\tau (\kern0.7pt y_{vo})/\rho }$ (see text). The abscissa $\tilde {y}^+$ represents the wall-normal coordinates in wall units determined using either the inner (left set of profiles for $y^+_{{in}} <5$) or the outer friction velocity (right set of profiles for $y^+_{{out}} > 50$) and considering either the wall or the virtual origin as the zero of the wall-normal axis: $\tilde {y}^+=yu_{\tau,in}/\nu$ and $\tilde {y}^+=(\kern0.7pt y-y_{vo})u_{\tau,out}/\nu$. The solid line without symbols refers to the profile of a canonical smooth wall channel flow at $Re_{\tau }=950$ (Hoyas & Jiménez Reference Hoyas and Jiménez2008). The sub-plot in panel (a) is a magnification of the profiles in the logarithmic region. (b) Roughness function (Hama Reference Hama1954) plotted as a function of $\lambda$ for all considered canopies. The solid line represents a polynomial fit for cases obtained varying $\Delta S$ (current investigation), while the dashed line represents a polynomial best fit for the cases considered by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). See table 1 for the symbols’ meaning.

Outside the canopy, the mean velocity profile is scaled with the outer friction velocity $u_{\tau,{out}}$ and the corresponding viscous length scale $\delta _{\nu,{out}} = u_{\tau,{out}}/\nu$. It follows a logarithmic distribution, and the presence of the canopy causes a shift in the logarithmic region (see inset in figure 11a). In this case, the logarithmic law for turbulent boundary layers over generic textured surfaces can be expressed as

(3.4)\begin{equation} U_{{out}}^+= \frac{1}{\kappa} \log(\kern0.7pt y_{{out}}^+) + 5.5 - \Delta U^+. \end{equation}

The first three terms in (3.4) represent the canonical log-law for a smooth wall, while the last term $\Delta U^+$ corresponds to the roughness function (Hama Reference Hama1954; Jiménez Reference Jiménez2004). The latter accounts for the momentum deficit or surplus due to the presence of textures on the wall and determines the wall offset in the logarithmic region. Examples of $\Delta U^+$ variations can be found in studies on permeable substrates (Rosti, Brandt & Pinelli Reference Rosti, Brandt and Pinelli2018; Gómez-de Segura & García-Mayoral Reference Gómez-de Segura and García-Mayoral2019) and rough wall boundary layers (Bechert, Bruse & Hage Reference Bechert, Bruse and Hage2000; Jiménez Reference Jiménez2004).

Figure 11(b) shows the roughness function as a function of the solidity $\lambda$ for all the canopies considered. All the canopy configurations considered in this work exhibit a drag increment, with the magnitude of the roughness function depending on $\lambda$. A polynomial fit is provided for the cases where $\Delta S$ is varied in the current investigation, while a dashed line represents the best fit for the cases considered by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). The symbols used in the figure are described in table 1.

Next, figure 12(a,b) present the mean filament spacing $\Delta S^+$ and the effective canopy height $k_{{eff}}^+=(h-y_{vo}) u_{\tau,{out}}/\nu$ as functions of $\lambda$ for all the considered canopies. The figure illustrates a distinct behaviour for solidity values obtained by varying the canopy height compared with those obtained by changing the mean spacing. The distributions of $\Delta S^+$ and $k_{{eff}}^+$ for the canopies considered by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) show a monotonous increase with $\lambda$, with a possible saturation of $\Delta S^+$ for very dense cases. However, when increasing the filament packing, a decrease in both $\Delta S^+$ and $k_{{eff}}^+$ is observed for increasing values of $\lambda$. The effective canopy height and mean spacing for the densest case (i.e. case HD) are found to be $\Delta S^+\approx 50$ and $k_{{eff}}^+\approx 20$, respectively, which correspond to the conditions on the protrusion height for the breakdown of the drag reduction regime in walls covered with riblets (Bechert et al. Reference Bechert, Bruse and Hage2000; Garcia-Mayoral & Jimenez Reference Garcia-Mayoral and Jimenez2011). A decreasing effective height and drag in the dense regime (i.e. when $\lambda >0.15$) can be related to the departure from a fully rough regime and the onset of a transitionally rough regime as previously reported by Jiménez (Reference Jiménez2004) and by Thakkar et al. (Reference Thakkar, Busse and Sandham2018). The transitional regime is characterised by a complex dynamical interplay between the drag offered by the roughness elements and the viscous drag (Flack Reference Flack2018).

Figure 12. (a) Mean filament spacing and (b) effective canopy height expressed in outer wall units ($u_{\tau,out}$), as a function of $\lambda$ for all considered canopies. See table 1 for symbols.

The equivalent sand roughness $k_s$ is another parameter that can be used to characterise the effect of non-smooth walls (Schlichting Reference Schlichting1937; Bradshaw Reference Bradshaw2000; Jiménez Reference Jiménez2004). Here, $k_s$ was introduced as a notional measure of the height of closely packed sand grains that would produce the same frictional drag in the fully rough regime (when $k_s^+ \geq 70$) (Nikuradse et al. Reference Nikuradse1950). Based on the collapse of the roughness function in the fully rough asymptotic condition (Flack & Schultz Reference Flack and Schultz2014), $k_s$ can be defined via a modified law of the wall (Jiménez Reference Jiménez2004) as

(3.5)\begin{equation} U_{out}^+=\frac{1}{\kappa}\log\left(\frac{y-y_{vo}}{k_s}\right)+8.5. \end{equation}

The ratio of $k_s/k$ (where $k$ is the roughness element height) and the roughness solidity are strongly related as shown in figure 13(a), where the non-dimensional form of the equivalent sand roughness ($k_s/k$) is plotted as a function of the effective solidity $\lambda _{{eff}}=dk/\Delta S^2$ for all the canopy configurations considered here and by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). In the graphs, the value of $k_s/k$ has been rescaled with a filament drag coefficient $C_d$ (Jiménez Reference Jiménez2004) computed at the mid-location of the filaments (the drag at any Lagrangian node along each filament can be readily obtained within the framework of the immersed boundary method we use (Pinelli et al. Reference Pinelli, Naqavi, Piomelli and Favier2010) as the force per unit volume necessary to establish the zero-velocity condition). The figure indicates that all the considered canopies belong to the sparse k-type regime, with the canopy-dense cases (HD, DE, MDM and MDE) in a transitional condition. Figure 13(b) shows the relationship between $k_s/\tilde {k}$ and $\lambda$ (i.e. the $h$ based canopy solidity). In the sparse k-type regime, $k_s/\tilde {k}$ increases linearly with $\lambda$, reaching a maximum at the transition threshold ($\lambda =0.15$). In denser regimes, the value of $k_s/\tilde {k}$ decreases with increasing $\lambda$ due to the sheltering effect of the filaments.

Figure 13. Non-dimensional form of equivalent sand roughness ($k_s/k$) as a function of the (a) effective solidity ($\lambda _{eff}$) and (b) solidity ($\lambda$). The value of $k_s/k$ has been rescaled with a filament drag coefficient $C_d$ (Jiménez Reference Jiménez2004) computed at the mid-location of the filaments, where the localised flow is unaffected by the fixed and free ends. The dashed lines are from Jiménez (Reference Jiménez2004). See table 1 for symbols. Other symbols represents non-circular roughness elements: $*$, spanwise fences (Schlichting Reference Schlichting1937); ‘x’, spanwise fences (Webb, Eckert & Goldstein Reference Webb, Eckert and Goldstein1971); ${\unicode{x2721}}$, spanwise cylinders (Tani Reference Tani1987).

While the concept of equivalent sand roughness can be helpful in establishing a connection between the outer behaviour of canopy flow and traditional experiments conducted on rough walls, it is evident from our findings that solely examining the external canopy flow is inadequate for fully characterising the overall impact of the canopy on the flow development above the canopy tip. In fact, the flow behaviour above the canopy exhibits similarities to various roughness morphologies, which are related to the solidity parameter.

The diagonal Reynolds stresses are analysed and presented in figure 14, where they are scaled with the outer friction velocity. The profiles collapse in the region outside the canopy, while the inner layer exhibits a substantial influence of the in-plane solidity in all three components (streamwise, wall-normal and spanwise). The outer region profiles collapse similarly to smooth wall turbulent channel flows, confirming outer layer similarity. The spanwise component distributions show an increase in the magnitude of the internal peak moving towards the canopy bed as the canopy solidity increases. The streamwise component variations exhibit a non-monotonic behaviour, with a decrease in magnitude for sparser regimes and an abrupt increase for denser regimes. The increase is accompanied by a movement of the peak towards the canopy bed, indicating the presence of another mechanism associated with the onset of the dense regime.

Figure 14. Diagonal Reynolds stresses plotted as a function of the wall-normal coordinate in external units ($H$). The stresses are scaled with the outer friction velocity computed at the virtual origin (see (3.2)). The panels represent (a) the streamwise, (b) the wall-normal and (c) the spanwise components. The solid line without symbols refers to the profile of a canonical smooth wall channel flow at $Re_{\tau }=950$ (Hoyas & Jiménez Reference Hoyas and Jiménez2008). See table 1 for symbols.

The diagonal Reynolds stress profiles are further analysed using a local velocity scale based on the total stress:

(3.6)\begin{equation} u_{\tau,l}=\sqrt{\frac{\nu\partial\langle U\rangle / \partial y -\langle u'v'\rangle}{1-y/H}}. \end{equation}

This is the same scaling proposed by Sharma & García-Mayoral (Reference Sharma and García-Mayoral2020a) in the context of canopy flows, and by Jiménez (Reference Jiménez2013) and Tuerke & Jiménez (Reference Tuerke and Jiménez2013) when dealing with manipulated wall-bounded turbulent flows. The profiles scaled with this local friction velocity (presented in figure 15) show a collapse across all cases, resembling smooth-wall-like behaviour, particularly in the canopy region. They exhibit an almost bimodal distribution, with two peaks located in the inner and outer regions. The magnitude of the inner peak increases with $\lambda$ for all three components. A detailed examination of the wall-normal component reveals that while the peak magnitude increases with $\lambda$ in the inner layer, it remains unaffected by the stem number density. The spanwise component profiles show an increase in the magnitude of the internal peak with increasing $\lambda$. The streamwise component variations exhibit an interesting behaviour, with a decrease in magnitude for sparser regimes and an increase for denser regimes, indicating the influence of the canopy density.

Figure 15. Diagonal Reynolds stresses, scaled with the local friction velocity, plotted as a function of the wall-normal coordinate in external units ($H$). The insets on the top right corner of panels (a–c) show the distributions within by the canopy layer. Panels (a–c) represent the streamwise, the wall-normal and spanwise components, respectively. Panel (d) shows the magnitude of the internal peak of the streamwise fluctuations as a function of $\Delta S/h$. See table 1 for symbol meanings.

3.4. Structure of the flows

We now turn our attention to the inception and arrangement of coherent structures within both the intra-canopy and outer regions, with a particular focus on the role of the canopy's in-plane solidity. We begin by examining the time-averaged flow field using a spatial ensemble averaging technique. Subsequently, we delve into the fluctuating flow field through triple decomposition and investigate the energy content of velocity fluctuations using spectral decomposition and 3-D autocorrelation functions.

The time averaged flow fields are obtained following a three-stages procedure (Monti et al. Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). First, we compute a time-averaged field encompassing the entire canopy. This mean field is not particularly meaningful as it contains local variations resulting from the random distribution of filaments across the tiles. To mitigate this effect, we proceed to the second stage, where we introduce a virtual cuboid with a base measuring $2\Delta S \times 2 \Delta S$ and a height equal to $h$. In this step, we adjust the velocity fields over each tile volume by shifting them to align their respective stems with the centre of the cuboid base. All the resulting translated fields, occupying the cuboid, are then ensemble-averaged to generate an intermediate mean field. Finally, in the third stage, we obtain the double-averaged field over a $\Delta S \times \Delta S$ tile, akin to what would correspond to a uniform filament distribution, by averaging across the $x$ and $z$ directions of the cuboid. This is achieved by leveraging the averaged velocity field's periodicity and parity.

Figure 16 presents the mean flow fields obtained using this procedure, displaying the mean streamlines and iso-contours of the wall-normal velocity. For visual clarity, the filaments located at $(x/\Delta S, z/\Delta S) = (-1, 0)$ and $(0, -1)$ in the HD configuration have been removed. In the sparser regime (figure 16a,b), the plane cut at the canopy tip reveals a mean ejection and sweep on the frontal and rear sides of the filaments. However, deeper within the canopy, we observe an opposite behaviour with the mean ejection and sweep occurring on the rear and frontal sides of the filaments at $y/h \approx 0.8$, indicating that the outer flow penetrates the canopy by sweeping along the rear side of each filament. Closer to the wall, the mean flow pattern becomes independent of the wall-normal coordinate.

Figure 16. Time and tile-ensemble averaged flow fields (see text for description of the averaging procedure) in the intra-canopy region corresponding to the five canopy configurations. Panels (a–e) correspond to the SP, TR, MD, DE and HD cases, respectively. The iso-contours represented in the cross-planes and wall-parallel planes, correspond to the wall-normal velocity, with blue and red regions representing negative and positive values within the range $[-0.3, 0.8]/U_b$. The wall-parallel planes correspond to the wall-normal locations of $5\,\%$, $45\,\%$ and $95\,\%$ of the canopy height. The streamlines have been computed using the tile-ensemble averaged flow field.

Furthermore, it is observed that the footprint of the mean sweeps observed near the canopy tip progressively diminishes as the canopy density increases (as seen in figure 16b–d). The sparser configurations, especially those in the mildly dense regime (i.e. MD and DE), exhibit very weak sweeps on the rear sides of the filaments, suggesting that the mean flow field is nearly unaffected by changes in the wall-normal coordinate.

An interesting observation arises from the mean flow analysis of the HD configuration (see figure 16c). It reveals the presence of a vortex between adjacent filaments, suggesting a decoupling phenomenon within denser canopy regions. This decoupling inhibits direct kinematic communication between the inner and outer flows, resulting in the emergence of a vortex with a size determined by filament spacing and local Reynolds number. This phenomenon bears resemblance to observations made in previous studies on flows over k-type rough surfaces (Perry et al. Reference Perry, Schofield and Joubert1969; Jiménez Reference Jiménez2004; Leonardi et al. Reference Leonardi, Orlandi and Antonia2007), further emphasising the similarity between the present HD configuration and these rough surface conditions.

In the mildly dense regime (MD and DE configurations), the mean flows resemble those of the fully dense configuration (MDE) studied by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). Partial decoupling is observed between the internal and external flows, limiting the formation of a stable vortex. The emergence of a mean vortex in the HD configuration can be conceptualised as a cavity-like flow, following principles proposed by Prasad & Koseff (Reference Prasad and Koseff1989) and Hou et al. (Reference Hou, Zou, Chen, Doolen and Cogley1995). This similarity becomes evident by looking at figure 17(a) that reports the rescaled, averaged streamwise and wall-normal velocities using the internal friction velocity ($u_{\tau,in}=(\tau _w/\rho )^{1/2}$) for the densest case within the void between successive filaments.

Figure 17. Time averaged velocity distributions within the canopy void centreline for the HD: dashed line is the streamwise component; solid line is the wall-normal one.

Cavity flows are influenced by non-dimensional parameters such as the geometric aspect ratio ($\Delta S/h$) and Reynolds number. The high shear at the canopy tip drives this cavity-like flow, with the Reynolds number linked to the local filament Reynolds number shown in figure 17(b) for various solidity values. It is noticed that the filament Reynolds number decreases with increasing stem density.

The velocity distribution of the streamwise velocity components shown in figure 17(a) also highlights the presence of a sheltering effect from the upstream filament. This effect (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004) is confirmed by the mean streamlines on the wall-parallel plane near the canopy tip (figure 16). In the sparser configurations, the streamlines show deflections around the filaments rejoining downstream. Differently, in denser canopies, the mean streamlines on the rear side of filaments do not reconnect, indicating a stronger sheltering effect. This is particularly evident in the HD case. Increased sheltering leads to reduced $\Delta U^+$ (figure 11b), as observed in studies on rough surfaces with densely packed elements (Jiménez Reference Jiménez2004).

We now shift our attention to investigating a potential instability within the canopy, as suggested by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). To explore this phenomenon, we examine figure 18, which presents instantaneous snapshots of velocity fluctuations in wall-parallel planes at locations corresponding to internal inflection points. In sparser configurations (SP and TR), the velocity fluctuations exhibit coherent regions resembling external streamwise-oriented vortices, along with smaller-scale fluctuations induced by fluid meandering around the filaments. However, in denser canopies, more pronounced modulation in the internal flow field becomes apparent. Panels ($\,j$) and ($m$) highlight strong spanwise-oriented structures resembling streamwise-moving fronts. Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) proposed that these spanwise coherent motions, reminiscent of Kelvin–Helmholtz-like instabilities, arise from high wall-normal permeability that promotes the directional penetration of external structures. Within the canopy, these penetrations generate streamwise and spanwise modulations due to bed impermeability and solenoidal conditions.

Figure 18. Instantaneous snapshots of the velocity fluctuations in the wall-parallel plane located at the inner inflection point. Panels (a,d,g,j,m), (b,e,h,k,n) and (c,f,i,l,o) represent the three velocity components $u'$, $v'$ and $w'$, respectively. The panels in the rows (a–c) to (m–o) display the cases ordered by increasing solidity. The iso-contours have been scaled with the local friction velocity. The colour codes, ranging from red to blue, represent positive to negative values. The velocity ranges are $u'^+ \in [-3,3]$, $v'^+ \in [-2,2]$ and $w'^+ \in [-3,3]$.

The variation and organisation of coherent motions within the canopy layer can be quantified through higher-order statistics of the full fluctuating field, specifically, the two-dimensional two-point auto-correlation function in the streamwise ($\Delta x$) and spanwise ($\Delta z$) directions,

(3.7)\begin{equation} R_{u'_iu'_i}=\frac{u_i'(0,y,0)u_i'(\Delta x,y,\Delta z)}{u_i'^{2}(0,y,0)}, \quad i=1,2, \end{equation}

yields the structure of the velocity fields associated with various regimes. Figure 19 illustrates the positive ($0.1$) and negative ($-0.1$) correlation coefficients, showing alterations in the red iso-surface related to the wall-normal velocity component correlation, except for the HD configuration, which exhibits a spiky structure. Conversely, the negatively correlated iso-surface (blue region) shrinks for larger $\lambda$ values until disappearing in the densest case.

Figure 19. Iso-surface of the two-dimensional two-point velocity auto-correlation function shown as a function of the wall-normal distance in the intra-canopy region. Panels (a,d,g,j,m), (b,e,h,k,n) and (c,f,i,l,o) correspond to the three velocity components, streamwise, wall-normal and spanwise, respectively. The panels in the rows (a–c) to (m–o) display the cases ordered by increasing solidity. The red and blue iso-surfaces represent the regions of positive and negative correlations. The green translucent slice shows the mean location of the inner inflection point while the yellow one indicates the position of the virtual origin.

The streamwise correlation function reveals that an increase in the in-plane solidity modifies the coherence length and direction, while the SP case shows long-scale correlation above and across the canopy. This behaviour links to the quasi-streamwise vortices that penetrate the canopy through sweep events.

As the canopy density increases, the penetration of external vortices into the canopy is governed by $\Delta S$. The HD configuration shows the shortest streamwise coherence length at the canopy top, where the vortices lose their coherence due to deformation. Furthermore, streamwise coherence decreases with proximity to the bottom wall, while spanwise coherence becomes more prominent.

The fluctuating flow field structure can be further explored by considering the spectral energy content of the fluctuating velocity components. The one-dimensional premultiplied spectra, rescaled with the local friction velocity, are shown in figures 20 and 21. They feature two peaks corresponding to the wavelengths $\lambda _x/H\simeq 1$ and $\lambda _z/H\simeq 1$, which are connected to the coherent structures populating the external flow.

Figure 20. Premultiplied spectra of the fluctuating velocity components represented as a function of the streamwise wavelength ($\lambda _x/H$) and wall-normal distance ($\kern0.7pt y/H$). The rows left to right refer to the three velocity components, with the first column (a,d,g,j,m), second column (b,e,h,k,n) and third column (c,f,i,l,o) representing the premultiplied spectra $\kappa _x\varPhi _{u'u'}/u_{\tau,l}^2$, $\kappa _x\varPhi _{v'v'}/u_{\tau,l}^2$ and $\kappa _x\varPhi _{w'w'}/u_{\tau,l}^2$, respectively. Solidity values increase moving from top to bottom. The yellow, red and green vertical lines correspond to the wavelengths matching the size of the diameter, the filament height and the spacing respectively. The cyan, yellow and red horizontal dashed lines represent the wall-normal location of the virtual origin, the location of the internal inflection point and the canopy tip.

Figure 21. Premultiplied spectra of the fluctuating velocity components $\kappa _z\varPhi _{u'u'}/u_{\tau,l}^2$, $\kappa _z\varPhi _{v'v'}/u_{\tau,l}^2$ and $\kappa _z\varPhi _{w'w'}/u_{\tau,l}^2$ represented as a function of the spanwise wavelength ($\lambda _z/H$) and wall-normal distance ($\kern0.7pt y/H$). For the panels’ ordering and meaning of the coloured lines, refer to the caption of figure 20.

The internal spectral peak ($\kern0.7pt y/H \leq 0.1$) points towards two scenarios arising in different canopy flow regimes. In the sparser regimes, the inner peaks do not correspond to the canopy geometrical parameters’ wavelengths. For larger $\lambda$, a peak emerges corresponding to the filament's diameter. This peak is most intense for the HD configuration, highlighting the wake's significance behind the filament in the dense regime.

In denser configurations, the internal peak of the spectral energy content splits into two maxima. The leftmost peak corresponds to the average filament spacing ($\lambda _x/H=\lambda _z/H=\Delta S$), while the rightmost mirrors the outer peak.

Finally, Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) associated the peaks of $u'$ and $w'$ with the energetic sweeps of the outer layer vortices into the canopy and subsequent deflection upon reaching the bottom wall, affected by the wall's impermeability and solenoidal constraints. The HD case exhibits a different behaviour due to the increased filament clustering affecting the wall-normal permeability, thus reducing the outer vortices’ role in determining the intra-canopy region's dynamics.

The internal structure of the intra-canopy layer is elucidated by considering the spectral decomposition of the fluctuating flow field. Figure 22 displays the two-dimensional premultiplied spectra in log coordinates, differentiated by colours. The red, blue and yellow regions denote the inner region, outer layer and overlap zone, respectively. Red and green lines illustrate wavelengths linked to the canopy height and mean filament spacing.

Figure 22. Iso-surfaces of two-dimensional premultiplied spectra of the fluctuating velocity components represented as a function of the wall-parallel wavelengths ($\lambda _x/H$ and $\lambda _z/H$) and the wall-normal distance ($\kern0.7pt y/H$). The first column (a,d,g,j,m), the second column (b,e,h,k,n) and third column (c,f,i,l,o) represent the $\kappa _x\kappa _z\varPhi _{u'u'}/u_{\tau,l}^2$, $\kappa _x\kappa _z\varPhi _{v'v'}/u_{\tau,l}^2$ and $\kappa _x\kappa _z\varPhi _{w'w'}/u_{\tau,l}^2$ premultiplied spectra, respectively. The rows from top to bottom correspond to decreasing in-plane solidity. See the text for the meaning of the lines and colours. Note that the isosurfaces are open-ended in the wall-normal direction because the region above the canopy is not reported here.

In sparser conditions (panels a–c), iso-surfaces span the canopy, exhibiting a wide range of wavelengths from $\lambda _x/H=\lambda _z/H=1$ to those tied to canopy geometry. Distinctive peaks within the spectra can be attributed to large-scale vortices infiltrating the canopy layer and smaller contributions from meandering filament motion.

With increased canopy density, the spectra's iso-surfaces reshape more intricately. As the solidity $\lambda$ rises, there is a more defined scale separation and energy content changes based on the wall-normal distance from the canopy bed. An elongated peak appears within a wavelength window of $\lambda _x/H,\lambda _z/H \in [h,\Delta S]$ across all components of the premultiplied spectra, associating with velocity fluctuations due to filament interaction.

Interestingly, the $v'$ component's spectral footprint displays two elongated peaks. Peaks extend further with increasing $\lambda$ values. A mildly elongated peak is noticed in the $u'$ component in the TR case (panel d), possibly due to transitional effects present at this $\lambda$ value.

In the HD case, a pronounced scale separation occurs in the intra-canopy region, visible in panels (m,n,o). The $u'$ component exhibits three isolated peaks at varying wavelengths, while the $v'$ and $w'$ components display peaks that correlate to large- and small-scale contributions. This suggests a correlation between the three velocity fluctuation components within the canopy. Peaks of the $v'$ component progressively decrease with an increasing $\lambda$.

One-dimensional premultiplied spectra do not clearly identify coherent internal motions. Large-scale contributions are confined to the canopy's upper region in the HD case, while small-scale contributions reside in the wavelength window $\lambda _x/H,\lambda _z/H \in [h,\Delta S]$.

Finally, our analysis returns to the external boundary layer above the canopy. Turbulence in this region is known to provoke large-scale spanwise coherent structures, identified either through the two-dimensional spectral content of wall-normal velocity fluctuations above the canopy tip, or the outer peak of the premultiplied cospectra of $u'$ and $v'$ components shown in figure 23. Figure 24 portrays two-dimensional premultiplied spectra at $y_{out}^+=60$ for all solidity values. Structures with strong coherence in the spanwise direction are indicated by a spanwise elongated peak within a narrow band of streamwise wavelengths.

Figure 23. Magnitude of premultiplied cospectra for the streamwise and the wall-normal velocity fluctuations as a function of the streamwise and spanwise wavelengths. Panels (a–e) represent the streamwise wavelength for increasing values of $\lambda$. Panels (f–j) represent the streamwise wavelength ordered with $\lambda$ as in the panels (a–e) cases. The contours are extracted within the range $[0,0.4]$ using a $0.02$ increment.

Figure 24. Two-dimensional premultiplied spectra of the wall-normal velocity fluctuations extracted at $y_{out}^+=60$, as a function of the streamwise ($\lambda _x/H$) and spanwise ($\lambda _z/H$) wavelength, scaled with the open channel height $H$. Panels (a–e) are ordered by increasing the solidity value (i.e. SP to H cases). The red contours highlight the spectra evaluated at $y^+=60$ for a smooth channel with $Re_\tau =710$.

The Kelvin–Helmholtz-type instability, analogous to a plane mixing layer, arises due to an inflected mean velocity profile. Many have used this analogy to study coherent structures at the canopy tip (Raupach et al. Reference Raupach, Finnigan and Brunet1996; Finnigan Reference Finnigan2000; Nepf Reference Nepf2012). Authors hypothesised that the spanwise roller wavelength ($\varLambda _x$) in a dense canopy, scaled with the shear length ($L_s$), falls within $7<\varLambda _x/L_s<10$. Raupach et al. (Reference Raupach, Finnigan and Brunet1996) proposed a precise estimate of $\varLambda _x \simeq 8.1L_s$.

However, the discrepancy between this estimate and figure 25(a), which presents wavelengths for all canopy configurations, suggests that the approximation holds mainly for tall, dense canopies. Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020) stated that canopy-interaction inconsistencies arise in sparser regimes, such as the MSP and MTR cases. Our simulations show that the wavelength for the sparsest configuration (SP case) aligns closely with the prediction of Raupach et al. (Reference Raupach, Finnigan and Brunet1996).

Figure 25. (a) Streamwise wavelength of the large spanwise coherent motions (made dimensionless with the canopy height $h$) above the canopy tip shown as a function of the shear length scale. The solid lines correspond to the current simulations, while the dashed line corresponds to the estimate proposed by Raupach et al. (Reference Raupach, Finnigan and Brunet1996). (b) Velocity difference within the canopy outer layer scaled by the bulk velocity, shown as a function of the solidity. Solid lines were obtained in the current investigation, and dashed lines were obtained from the study by Monti et al. (Reference Monti, Omidyeganeh, Eckhardt and Pinelli2020). See table 1 for symbols.

The inconsistency can be attributed to differing formulations. In the formulation by Raupach et al. (Reference Raupach, Finnigan and Brunet1996), the velocity profile reaches negligible value within a distance $L_s$ from the tip, which does not account for shorter canopies. To include sparser canopies, we have adjusted the calculation to $U_o=U_h-L_s (dU/{{d} y})_{y=h}$ using the mean velocity at the virtual origin ($U_h$ being the velocity at the canopy tip, i.e. at $y=h$).

Our findings, presented in figure 25(b), indicate that shear length scale's validity as vorticity thickness is constrained to canopies with large velocity differences, specifically those exceeding a threshold magnitude of $\Delta U_s/U_b \approx 0.26$.