1 Introduction
The present work studies flows over dense canopies of filaments of small size. Canopy flows are mainly studied in the context of natural vegetation canopies, but they also encompass engineering flows where the canopy parameters may be very different from those of natural canopies. Many of the key findings from natural canopy studies have been summarised in the reviews by Finnigan (Reference Finnigan2000), Belcher, Harman & Finnigan (Reference Belcher, Harman and Finnigan2012) and Nepf (Reference Nepf2012). In engineering applications, filament canopies can, for instance, be used to enhance heat transfer (Fazu & Schwerdtfeger Reference Fazu and Schwerdtfeger1989; Bejan & Morega Reference Bejan and Morega1993) and for energy harvesting (McGarry & Knight Reference McGarry and Knight2011; Elahi, Eugeni & Gaudenzi Reference Elahi, Eugeni and Gaudenzi2018). Depending on the geometry and spacing of their elements, canopies can be classified as sparse, dense or transitional (Nepf Reference Nepf2012). Dense canopies typically have small element spacings compared to the length scales in the overlying flow, and thus prevent turbulent eddies from penetrating efficiently within the canopy. Sparse canopies, on the other hand, have large element spacings and consequently, turbulent eddies are essentially able to penetrate the full height of the canopy (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004; Nepf Reference Nepf2012; Sharma & García-Mayoral Reference Sharma and García-Mayoral2018, Reference Sharma and García-Mayoral2020). Transitional, or intermediate, canopies would lie between these two regimes. In the present study, we assess how the flow within and above dense canopies is affected by canopy parameters, such as the element height and spacing. The canopies considered have spacings 
$s^{+}\approx O(10)$, which should be small enough to limit the penetration of the overlying turbulence within them. The frontal area density 
$\unicode[STIX]{x1D706}_{f}$ is also a commonly used measure to categorise canopies (Finnigan Reference Finnigan2000; Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004; Huang, Cassiani & Albertson Reference Huang, Cassiani and Albertson2009; Nepf Reference Nepf2012). Canopies with 
$\unicode[STIX]{x1D706}\gg 0.1$ are classified as dense, with 
$\unicode[STIX]{x1D706}\approx 0.1$ as transitional and with 
$\unicode[STIX]{x1D706}\ll 0.1$ as sparse. However, in addition to 
$\unicode[STIX]{x1D706}_{f}$, the length scales of the overlying turbulence should also be considered when determining the canopy regime. A given canopy geometry with a fixed 
$\unicode[STIX]{x1D706}_{f}$ may have element spacings much smaller than any overlying turbulent eddy at a particular Reynolds number, thereby not allowing turbulence to penetrate within the canopy. As the Reynolds number is increased, however, the size of these eddies will eventually become comparable to the element spacing, allowing turbulence to penetrate efficiently within the canopy. To assess this effect, we also study canopies with self-similar geometries, which have a fixed 
$\unicode[STIX]{x1D706}_{f}$, but different sizes in friction units.
We also place attention on the effect of canopy parameters on the Kelvin–Helmholtz-like, mixing-layer instability characteristic of dense canopy flows (Raupach, Finnigan & Brunet Reference Raupach, Finnigan and Brunet1996; Finnigan Reference Finnigan2000; Nepf Reference Nepf2012). This instability originates from the inflection point in the mean velocity profile at the canopy-tip plane (Raupach et al. Reference Raupach, Finnigan and Brunet1996). Kelvin–Helmholtz instabilities manifest as spanwise coherent rollers whose streamwise scale is determined by the shear-layer thickness (Michalke Reference Michalke1972; Brown & Roshko Reference Brown and Roshko1974). Ghisalberti & Nepf (Reference Ghisalberti and Nepf2004) noted that, while in free-shear flows the shear-layer thickness, and consequently the instability wavelength, continues to grow downstream, in fully developed canopy flows this thickness is constant and is set by the net canopy drag. Therefore, a fixed instability wavelength is generally associated with dense canopy flows. Several studies have shown that some aspects of this instability can be captured using a mean-flow linear stability analysis (Raupach et al. Reference Raupach, Finnigan and Brunet1996; White & Nepf Reference White and Nepf2007; Singh et al. Reference Singh, Bandi, Mahadevan and Mandre2016; Zampogna et al. Reference Zampogna, Pluvinage, Kourta and Bottaro2016; Luminari, Airiau & Bottaro Reference Luminari, Airiau and Bottaro2016). Some studies have also suggested that at the high Reynolds numbers of natural canopy flows these instabilities can be distorted by the ambient turbulence fluctuations and lose their spanwise-coherent nature (Finnigan, Shaw & Patton Reference Finnigan, Shaw and Patton2009; Bailey & Stoll Reference Bailey and Stoll2016). The importance of this instability decreases as the element spacing is increased, and sparse canopies do not exhibit a notable signature (Poggi et al. Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004; Huang et al. Reference Huang, Cassiani and Albertson2009; Pietri et al. Reference Pietri, Petroff, Amielh and Anselmet2009; Sharma & García-Mayoral Reference Sharma and García-Mayoral2020).
Based on the observations of previous studies, we would expect that the effect of increasing the canopy height for a fixed element spacing on the instability and the surrounding flow would eventually saturate. Ghisalberti (Reference Ghisalberti2009) and Nepf (Reference Nepf2012) proposed that the effective canopy height perceived by the overlying flow would be a function of the canopy shear-layer thickness as it determined the extent to which the Kelvin–Helmholtz-like instabilities penetrated within the canopy. In flows over arrays of cuboidal posts, Sadique et al. (Reference Sadique, Yang, Meneveau and Mittal2017) found that the mean-velocity profiles over them became independent of the element heights at large element aspect ratios. They concluded that the overlying flow only interacted with the region near the element tips, and that the height below this ‘active’ region was dormant, and did not have a significant effect on the overlying flow. For their geometries, Sadique et al. (Reference Sadique, Yang, Meneveau and Mittal2017) observed the height of this active region to be related to the element width. A similar observation was also made by MacDonald et al. (Reference MacDonald, Ooi, García-Mayoral, Hutchins and Chung2018), who performed direct numerical simulations (DNS) of flows over spanwise-aligned bars. They found that the gap between the bars was the relevant length scale for the overlying flow, and that increasing the height of the bars beyond a certain height-to-gap ratio did not affect the overlying flow, or cause an increase in the drag they produced.
In the present work, we conduct a systematic range of DNS changing the canopy height and spacing separately in order to study their individual effects on the surrounding turbulence and on the Kelvin–Helmholtz-like instability. We also consider canopy geometries with constant 
$\unicode[STIX]{x1D706}_{f}$ for which the height and spacing are changed simultaneously in a fixed proportion. The canopies consist of rigid, prismatic filaments with small element spacings and large height-to-spacing ratios. The element spacings considered, 
$s^{+}\approx 3{-}50$, are much smaller than those typical of most natural canopy flows and would, for instance, be representative of flows over engineered canopies such as those mentioned previously in this section. We also assess how models based on linear stability analysis capture some of the effects of the canopy parameters on the Kelvin–Helmholtz-like instability.
The paper is organised as follows. The numerical methods used for the simulations and the canopy parameters are discussed in § 2. The results from the DNS, detailing the effect of the canopy parameters on the overlying turbulence and the Kelvin–Helmholtz-like instabilities are discussed in § 3. The results from linear stability analysis and a model to capture the instabilities are presented in § 4. The conclusions are summarised in § 5.
2 Methodology
We conduct DNS of symmetric channels with rigid canopy elements on both walls. The streamwise, wall-normal and spanwise coordinates are 
$x$, 
$y$ and 
$z$, with the associated velocities 
$u$, 
$v$ and 
$w$, and 
$p$ is the kinematic pressure. The wall-normal origin, 
$y=0$, is defined at the tip plane of the canopies protruding from the bottom wall. The channel height, 
$2\unicode[STIX]{x1D6FF}$, is defined as the distance between the tip planes of the canopies on the top and bottom walls. The canopy elements, therefore, extend below 
$y=0$ and above 
$y=2\unicode[STIX]{x1D6FF}$ and have a height 
$h$. A schematic representation of the channel is portrayed in figure 1. The size of the domain is a standard 
$2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ in the streamwise direction and 
$\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ in the spanwise direction. We use the channel half-height as the length scale in outer units, which implies that 
$\unicode[STIX]{x1D6FF}=1$ in outer scaling. The domain-to-canopy height ratio for most cases considered here is 
$(\unicode[STIX]{x1D6FF}+h)/h\approx 3$. We will show in § 3 that the height of the roughness sublayer scales with the canopy spacing rather than their height, as in the configurations of Sadique et al. (Reference Sadique, Yang, Meneveau and Mittal2017) and MacDonald et al. (Reference MacDonald, Ooi, García-Mayoral, Hutchins and Chung2018), and that outer-layer similarity is recovered well below the channel half-height. The channel height to element spacing ratio for most canopies considered is 
$\unicode[STIX]{x1D6FF}/s\gtrsim 10$, and for the canopy with the largest element spacing is 
$\unicode[STIX]{x1D6FF}/s\approx 4$. The flow is incompressible and the density is always scaled with the fluid density, implying that 
$\unicode[STIX]{x1D70C}=1$. The simulations are run at a constant flow rate, with the viscosity, 
$\unicode[STIX]{x1D708}$, adjusted to obtain a friction Reynolds number 
$Re_{\unicode[STIX]{x1D70F}}=u_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D708}\approx 185$ for most of the cases, where 
$u_{\unicode[STIX]{x1D70F}}$ is the friction velocity calculated at the canopy tips. In order to ascertain the effects of the Reynolds number on the flow, a simulation at 
$Re_{\unicode[STIX]{x1D70F}}\approx 405$ was also conducted. The simulation parameters are given in table 1 for reference. Scaling with 
$u_{\unicode[STIX]{x1D70F}}$ and 
$\unicode[STIX]{x1D708}$ is referred to as in friction or wall units, and scaling with the channel bulk velocity, 
$U_{b}$, and 
$\unicode[STIX]{x1D6FF}$ is referred to as in outer units.
Schematic representation of the domain considered in the present study.
Simulation parameters. Here, 
$N_{x}$ and 
$N_{z}$ are the number of rows of canopy elements in the streamwise and spanwise directions, respectively. The number of points used to resolve each period of the canopy in the streamwise and spanwise directions are 
$n_{x}$ and 
$n_{z}$, respectively. Here, 
$u_{\unicode[STIX]{x1D70F}}$ is the friction velocity based on the shear at the canopy tips scaled with the channel bulk velocity, 
$Re_{\unicode[STIX]{x1D70F}}$ is the friction Reynolds number based on 
$u_{\unicode[STIX]{x1D70F}}$ and 
$\unicode[STIX]{x1D6FF}$. The canopy frontal area density, height, spacing and width are 
$\unicode[STIX]{x1D706}_{f}$, 
$h$, 
$s$ and 
$w$, respectively.
The numerical method used to solve the three-dimensional Navier–Stokes equations is adapted from Fairhall & García-Mayoral (Reference Fairhall and García-Mayoral2018). A Fourier spectral discretisation is used in the streamwise and spanwise directions. The wall-normal direction is discretised using a second-order centred difference scheme on a staggered grid. The grid in the wall-normal direction is stretched to give a resolution 
$\unicode[STIX]{x0394}y_{min}^{+}\approx 0.33$ at the canopy-tip plane, stretching to 
$\unicode[STIX]{x0394}y_{max}^{+}\approx 3.3$ at the channel centre. The grid within the canopies preserves the resolution of 
$\unicode[STIX]{x0394}y_{min}^{+}\approx 0.33$ near the canopy-tip plane, and for the tallest canopies considered stretches to 
$\unicode[STIX]{x0394}y_{max}^{+}\approx 4$ at the base of the canopy, where the flow is quiescent. The wall-normal grid distribution for a representative canopy simulation is provided in appendix A for reference. To resolve the element-induced flow while avoiding excessive computational costs, the domain is divided into three blocks in the wall-parallel directions (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011; Fairhall & García-Mayoral Reference Fairhall and García-Mayoral2018; Abderrahaman-Elena, Fairhall & García-Mayoral Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019). In the central block, the resolutions in the streamwise and spanwise directions are 
$\unicode[STIX]{x0394}x^{+}\approx 6$ and 
$\unicode[STIX]{x0394}z^{+}\approx 3$, respectively, sufficient to resolve the turbulent eddies. The blocks including the canopy elements and the roughness sublayer have a finer resolution than the central block. In the fine blocks, the limiting resolution is not the one required to resolve the turbulent scales, but that required to resolve the obstacles or the element-induced flow. The resolutions in these blocks are given in table 1. The height of the fine blocks is chosen such that the element-induced flow decays to zero well within the fine-block region, and this is verified a posteriori. The time advancement is carried out using a three-step Runge–Kutta method with a fractional step, pressure correction method to enforce continuity (Le & Moin Reference Le and Moin1991)
$$\begin{eqnarray}\displaystyle \left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}\right]\boldsymbol{u}_{k}^{n} & = & \displaystyle \boldsymbol{u}_{k-1}^{n}+\unicode[STIX]{x0394}t\left[\frac{\unicode[STIX]{x1D6FC}_{k}}{Re}\text{L}\boldsymbol{u}_{k-1}^{n}-\unicode[STIX]{x1D6FE}_{k}\text{N}(\boldsymbol{u}_{k-1}^{n})\right.\nonumber\\ \displaystyle & & \displaystyle -\left.\unicode[STIX]{x1D701}_{k}\text{N}(\boldsymbol{u}_{k-2}^{n})-(\unicode[STIX]{x1D6FC}_{k}+\unicode[STIX]{x1D6FD}_{k})\text{G}(p_{k}^{n})\right],k\in [1,3],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{DG}(\unicode[STIX]{x1D719}_{k}^{n})=\frac{1}{(\unicode[STIX]{x1D6FC}_{k}+\unicode[STIX]{x1D6FD}_{k})\unicode[STIX]{x0394}t}\text{D}(\boldsymbol{u}_{k}^{n}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{k+1}^{n}=\boldsymbol{u}_{k}^{n}-(\unicode[STIX]{x1D6FC}_{k}+\unicode[STIX]{x1D6FD}_{k})\unicode[STIX]{x0394}tG(\unicode[STIX]{x1D719}_{k}^{n}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p_{k+1}^{n}=p_{k}^{n}+\unicode[STIX]{x1D719}_{k}^{n}, & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D644}$ is the identity matrix and L, D and G are the Laplacian, divergence and gradient operators, respectively. Here, N is the advective term which is dealiased using the 
$2/3$-rule (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2012). The Runge–Kutta coefficients, 
$\unicode[STIX]{x1D6FC}_{k}$, 
$\unicode[STIX]{x1D6FD}_{k}$, 
$\unicode[STIX]{x1D6FE}_{k}$ and 
$\unicode[STIX]{x1D701}_{k}$, for each substep, 
$k$, are taken from Le & Moin (Reference Le and Moin1991). The time step is 
$\unicode[STIX]{x0394}t$.
The canopy elements are represented using an immersed-boundary method adapted from García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011). Further details about the immersed-boundary method and validation studies are provided in appendix A. The parameters of the different simulations conducted are summarised in table 1. The simulation denoted by ‘SC’ is of a turbulent channel flow with smooth walls. The canopy-flow simulations are divided into three groups. The canopy elements studied in each group are prismatic, with a square top-view cross-section, and their arrangement is illustrated in figure 2. The first group, denoted by the prefix ‘S’, consists of canopies with a fixed height, 
$h^{+}\approx 96$, and element spacings ranging from 
$s^{+}\approx 10$ to 48. The second group, marked by the prefix ‘H’, consists of canopies with a fixed element spacing, 
$s^{+}\approx 16$, and element heights ranging from 
$h^{+}\approx 16$ to 128. The element width-to-spacing ratio for the canopies of S and H is 
$w/s=1/2$. The final group, denoted by the prefix ‘G’, consists of self-similar elements with a fixed height-to-spacing ratio 
$h/s\approx 4$, and 
$w/s=2/9$. The heights for the canopies of G range from 
$h^{+}\approx 10$ to 100, with the element spacings varying in proportion to their height. These canopies have a constant 
$\unicode[STIX]{x1D706}_{f}=0.85$ and are used to study the effect of changing the canopy size for a fixed geometry. Two additional simulations, 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$, are conducted to check the dependence of the results on the friction Reynolds number. The canopy geometries for both these simulations have 
$s^{+}\approx 16$, 
$h^{+}\approx 32$ and 
$w/s=1/2$, with friction Reynolds numbers 
$Re_{\unicode[STIX]{x1D70F}}\approx 180$ and 400. We also conducted several simulations to assess whether the wall-parallel resolutions used in the simulations are sufficient to resolve the element-induced flow. The simulation S24 was run at resolutions of 12, 24 and 36 points per element spacing, and G100 at 9, 18 and 27 points per element spacing. Different resolution sets are used for the geometries of S and G as they have different element width-to-spacing ratios. The root mean square (r.m.s.) velocity fluctuations obtained from these simulations are portrayed in appendix A. The simulation results are grid independent at a resolution of 24 points per element spacing for the geometry of case S24, and 18 points per spacing for that of case G100. The simulations with 9 and 12 points per spacing tend to under-predict the fluctuations within the canopies, with the maximum deviation observed in the wall-normal fluctuations of 20 % within the canopy. This discrepancy reduces to 4 % outside the canopy. These resolutions are only used for the densest canopy cases, where the fluctuations within the canopy are already very small, and therefore, higher resolution simulations would not change the trends observed. The higher Reynolds number simulation, case 
$\text{H}32_{400}$, is also simulated using 12 points per spacing. For this simulation, using a higher resolution would be computationally restrictive. Note that the same resolutions are used for cases 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$ to avoid grid related discrepancies in the comparison of their results.
Schematic of the canopy layouts considered in the present study. The canopies are characterised by their element height, 
$h$, the element width, 
$w$, and the element spacing, 
$s$. Note that the elements have a square top-view cross-section.
2.1 Reynolds number effect
Root mean square velocity fluctuations and Reynolds shear stresses for cases 
$\text{H}32_{180}$ in red and 
$\text{H}32_{400}$ in blue. The black lines represent the corresponding smooth-wall cases. The data for the smooth-wall simulations at 
$Re_{\unicode[STIX]{x1D70F}}\approx 400$ are taken from Moser, Kim & Mansour (Reference Moser, Kim and Mansour1999).
To analyse the influence of the Reynolds number in our subsequent DNS, we compare the results of cases 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$, which have the same canopy height and spacings in friction units, but different friction Reynolds numbers. The velocity fluctuations and the Reynolds shear stresses within the canopy, and above it up to a height of 
$y^{+}\approx 10$, of these simulations essentially collapse, as shown in figure 3. This suggests that the flow in the region near the canopy-tip plane scales in friction units, similar to the near-wall region in smooth-wall flows (Moser et al. Reference Moser, Kim and Mansour1999). Scaling in friction units over conventional rough surfaces has also been noted by Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015). Beyond 
$y^{+}\gtrsim 10$, we observe that the magnitude of the peaks in the fluctuations and the Reynolds shear stresses are larger for case 
$\text{H}32_{400}$ compared to case 
$\text{H}32_{180}$. The increase in magnitude of the near-wall peaks in the velocity fluctuations at friction Reynolds numbers larger than 
$Re_{\unicode[STIX]{x1D70F}}\approx 180$ is consistent with that observed in smooth-wall flows (Moser et al. Reference Moser, Kim and Mansour1999; Sillero, Jiménez & Moser Reference Sillero, Jiménez and Moser2013), also included in figure 3 for reference. Further away from the canopy tips, at 
$y^{+}>50$, the r.m.s. velocity fluctuations from the canopy simulations coincide with those from the smooth-wall simulations at their corresponding Reynolds numbers, which indicates the recovery of outer-layer similarity. In addition to the r.m.s. fluctuations being similar for these simulations, the distribution of energy in different scales is also similar. This is illustrated by the premultiplied spectral energy densities at 
$y^{+}\approx 15$, portrayed in figure 4. This height roughly corresponds to the location where the magnitude of the fluctuations peaks in smooth-wall flows (Jiménez & Pinelli Reference Jiménez and Pinelli1999). The results of 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$ suggest that the effect of the canopy scales in friction units, and therefore the results presented in the following sections for flows at 
$Re_{\unicode[STIX]{x1D70F}}\approx 180$ should also be relevant for higher Reynolds number flows.
Premultiplied spectral energy densities for cases 
$\text{H}32_{180}$ (line contours) and 
$\text{H}32_{400}$ (shaded contours), normalised by the respective r.m.s. values, at a height 
$y^{+}\approx 15$. Contours from (a–d) are in increments of 0.075, 0.06, 0.07 and 0.1, respectively.
Instantaneous realisations of the wall-normal velocity at 
$y\approx 0.1s$, normalised by 
$u_{\unicode[STIX]{x1D70F}}$. From top to bottom, (a,c,e,g,i) represent cases H16 to H128; and (b,d,f,h,j), cases G10 to G100. The insets in (b) and (d) provide a magnified view of the region in the bottom left corner of these panels, marked with a black rectangle. The clearest and darkest colours represent intensity 
$\pm 0.4$ in (a,c,e,g,i) and, from top to bottom, 
$\pm (0.2,0.4,0.8,0.8,1.0)$ in (b,d,f,h,j).
Root mean square velocity fluctuations of the element-induced flow. The lines from red to blue, indicated by the direction of the arrows, represent (a,d,g) cases S10 to S48; (b,e,h) cases H16 to H128; and (c,f,i) cases G10 to G100.
3 Effect of canopy parameters on the surrounding turbulence
In this section, we discuss the results obtained from the DNS, aiming to characterise the effect of the canopy parameters on both the Kelvin–Helmholtz-like instability and the fluctuating flow within and above the canopies.
3.1 Height of the roughness sublayer
Before discussing the effect of the canopy on the overlying flow, we first define the extent of the region affected by the canopy, that is, the height of the roughness sublayer. Over conventional rough surfaces, with heights comparable to or smaller than their spacings, the height of the roughness sublayer is generally observed to be a function of the roughness height (Raupach, Antonia & Rajagopalan Reference Raupach, Antonia and Rajagopalan1991; Flack, Schultz & Connelly Reference Flack, Schultz and Connelly2007; Abderrahaman-Elena et al. Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019). Jiménez (Reference Jiménez2004) reviewed the effect of various roughness geometries on turbulent flows and noted that, in flows over closely packed spanwise aligned grooves, the flow within each groove would be isolated from the overlying flow due to the ‘sheltering’ effect of the preceding obstacle. The overlying flow in this case would not interact with the full height of the groove. This sheltering effect was also noted by Sadique et al. (Reference Sadique, Yang, Meneveau and Mittal2017) for high-aspect-ratio prismatic roughness and by MacDonald et al. (Reference MacDonald, Ooi, García-Mayoral, Hutchins and Chung2018) for spanwise aligned grooves with large spacings, and has been used to model cuboidal roughness by Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016). As the element spacings of the canopies studied here are small, this sheltering effect should result in the overlying flow only interacting with the region near the canopy-tip plane. In order to determine the height of this region, we examine the element-coherent flow induced by the canopy elements. The footprint of the element-induced flow can be observed in the instantaneous realisations of the velocity above the canopy-tip plane, portrayed for the canopies of families H and G in figure 5. We isolate the element-induced flow using the standard triple decomposition of Reynolds & Hussain (Reference Reynolds and Hussain1972)
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}=\boldsymbol{U}+\boldsymbol{u}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\prime }=\widetilde{\boldsymbol{u}}+\boldsymbol{u}^{\prime \prime }, & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{u}$ is the full velocity, 
$\boldsymbol{U}$ is the mean velocity obtained by averaging the flow in time and space, and 
$\boldsymbol{u}^{\prime }$ is the full space- and time-fluctuating signal. The latter is decomposed into the element-induced, dispersive velocity, 
$\widetilde{\boldsymbol{u}}$, which is obtained by ensemble-averaging the flow in time alone, and the element-incoherent fluctuating velocity 
$\boldsymbol{u}^{\prime \prime }$, which includes the contributions from the background turbulence and the Kelvin–Helmholtz-like instability. The r.m.s. fluctuations of 
$\widetilde{\boldsymbol{u}}$, therefore, result from fluctuations in space alone.
We observe that the element-induced fluctuations, for all the canopies studied here, decay exponentially above the canopy-tip plane, and become negligible at a height of one element spacing above regardless of the canopy depth, as shown in figure 6. This suggests that the height of influence of the element-induced flow is determined by the spacing between the elements rather than their height.
Even though the element-induced fluctuations only extend to one element spacing above the canopy-tip plane, their influence on the background turbulence extends to a height of approximately 2–3 element spacings, as can be observed in figure 7. Abderrahaman-Elena et al. (Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019) observed a similar effect over conventional cubical roughness, where the element-induced fluctuations only extended to 
$y\approx h$, but the effect of the roughness on the overlying flow extended to 
$y\approx 3h$ above them. At heights of 
$y/s>2{-}3$ above the canopy-tip plane, the full r.m.s. velocity fluctuations collapse with those of smooth-wall turbulence, as shown in figure 7, which is indicative of the recovery of outer-layer similarity. This is verified by a comparison of the premultiplied spectral energy densities of the canopy and smooth-wall cases in figure 8, which shows that the energy densities of the canopies of family S collapse with those of the smooth-wall case for 
$y^{+}\gtrsim 90$. This corresponds to a height of approximately 
$2s$ for case S48, the canopy with the largest spacing. Although not shown, the premultiplied spectral energy densities of the canopies of families H and G collapse with the smooth-wall spectra for 
$y/s\gtrsim 3$ as well. Previous canopy studies have proposed that the influence of the canopy elements on the flow above them is set by the wall-normal extent of the Kelvin–Helmholtz-like instability, which is determined by the canopy shear-layer thickness (Ghisalberti Reference Ghisalberti2009; Ghisalberti & Nepf Reference Ghisalberti and Nepf2009; Nepf Reference Nepf2012). It will be demonstrated in § 4 that the shear-layer thickness of the present canopies also depends on the element spacing.
Root mean square velocity fluctuations within and above the canopies. The lines from red to blue, indicated by the direction of the arrows, represent (a,d,g) cases S10 to S48; (b,e,h) cases H16 to H128; and (c,f,i) cases G10 to G100. The black lines represent the smooth-wall case, SC.
Premultiplied spectral energy densities at 
$y^{+}\approx 90$, with line contours from red to blue representing cases S10 to S48, normalised by their respective 
$u_{\unicode[STIX]{x1D70F}}$. The filled contours represent the smooth-wall case, SC. The contours in (a–d) are in increments of 0.11, 0.04, 0.06 and 0.04, respectively.
3.2 Effect of element height and spacing
In this section, we discuss the effect of the element height and spacing on the element-induced and on the full velocity fluctuations, both within and above the canopy. We observe that the element-induced fluctuations are largest near the canopy-tip plane and decay below it, as shown in figure 6, because they are obstructed by the canopy elements. This effect is more intense for smaller element spacings, and eventually results in the fluctuations vanishing completely well above the canopy floor for the simulations with the smallest spacings in families S and G. For the canopies of family H, which have a constant element spacing, the change in element height does not have a noticeable effect on the element-induced fluctuations, as shown in figure 6(b,e,h). For the canopies of families S and G, the intensity of the element-induced velocity fluctuations within the canopy increases with element spacing, when scaled using either the friction velocity or the channel bulk velocity. These results suggest that, for canopy elements with a given width, the magnitude of the element-induced fluctuations is governed mainly by the element spacing.
As discussed in § 3.1, in canopies with very small element spacing the height of the roughness sublayer is small, and we would expect such canopies not to disrupt the overlying turbulence significantly, regardless of their depth. In the literature on conventional roughness, small roughness elements that have a negligible effect on the overlying turbulent flow are termed ‘hydraulically smooth’, as the flow over them remains essentially smooth-wall like (Nikuradse Reference Nikuradse1933; Raupach et al. Reference Raupach, Antonia and Rajagopalan1991). Roughness elements with a characteristic size of a few wall units, 
$h^{+}\lesssim 5$, typically fall into the hydraulically smooth category (Raupach et al. Reference Raupach, Antonia and Rajagopalan1991; Jiménez Reference Jiménez2004; Flack et al. Reference Flack, Schultz and Connelly2007). Of the canopies studied here, we observe that the overlying flow for canopy G10, which has an element spacing of 
$s^{+}\approx 2.6$, is essentially smooth-wall-like above the canopy-tip plane. This is evidenced by the collapse of the r.m.s. velocity fluctuations, Reynolds shear stresses, and the mean velocity profile of this case with those of the smooth-wall case, as shown in figures 7(c,f,i) and 9(c,f). In addition, the magnitude of the velocity fluctuations below the canopy-tip plane is negligible. This suggests that the overlying turbulent flow essentially perceives the canopy-tip plane as an impermeable wall, and has little or no interaction with the canopy region below this plane.
For canopies with larger element spacings, we begin to observe deviations from smooth-wall-like behaviour in the overlying flow. Above the canopy-tip plane, an increase in the element spacing causes a reduction in the intensity of the streamwise velocity fluctuations and an increase in the intensity of the wall-normal and spanwise ones, as can be observed in figure 7 for the canopies of families S and G. For canopies with large element spacings, such as those of S48 and G100, the peak in 
$u^{\prime }$ typical of smooth-wall flows, is significantly reduced. These changes in the velocity fluctuations are accompanied by a reduction in the streamwise coherence in the flow with increasing element spacing as can be observed in the instantaneous realisations of the wall-normal velocity for the canopies of family G, portrayed in figure 5. Near-wall turbulence over smooth walls is characterised by streaks and quasi-streamwise vortices, which are predominantly streamwise-coherent (Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967; Jiménez & Pinelli Reference Jiménez and Pinelli1999). The decrease in 
$u^{\prime }$ and increase in 
$v^{\prime }$, 
$w^{\prime }$ above the canopy with increasing element size is also commonly reported over conventional rough surfaces (Ligrani & Moffat Reference Ligrani and Moffat1986; Orlandi & Leonardi Reference Orlandi and Leonardi2006). Several authors have attributed these changes in the velocity fluctuations and the loss of streamwise coherence to the roughness elements modifying the near-wall cycle and turbulence becoming more ‘isotropic’ (Jiménez Reference Jiménez2004; Flores & Jiménez Reference Flores and Jiménez2006; Flack et al. Reference Flack, Schultz and Connelly2007; Abderrahaman-Elena et al. Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019). We also observe an increase in the Reynolds shear stresses above the canopy tip plane with increasing element spacing, shown in figure 9(d,f), with an associated increase in the drag. The drag increase caused by rough surfaces is generally expressed in terms of the downward shift in the logarithmic region of the mean-velocity profile compared to that for a smooth wall (Hama Reference Hama1954). This shift can be observed for the canopies of families S and G in figure 9(a,c).
Profiles of the (a–c) streamwise mean velocity and (d–f) Reynolds shear stresses. The lines from red to blue, indicated by the direction of the arrows, represent (a,d) cases S10 to S48; (b,e) cases H16 to H128; and (c,f) cases G10 to G100. The black lines represent the smooth-wall case, SC.
Focusing now on the flow within the canopy, increasing the element spacing results in an increase in the magnitude of all the components of the full velocity fluctuations, as shown in figure 7, which is consistent with the observations of Green, Grace & Hutchings (Reference Green, Grace and Hutchings1995), Novak et al. (Reference Novak, Warland, Orchansky, Ketler and Green2000), Poggi et al. (Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004) and Pietri et al. (Reference Pietri, Petroff, Amielh and Anselmet2009). The wall-parallel velocity fluctuations, 
$u^{\prime }$ and 
$w^{\prime }$, decay rapidly below the canopy-tip plane, and their magnitude reaches a plateau in the core of the canopy, before dropping again near the canopy base to meet the no-slip condition. The abrupt changes in the velocity fluctuations near the element tips are typical of textures with perfectly flat and aligned tips and have also been observed over conventional cuboidal rough surfaces (Leonardi & Castro Reference Leonardi and Castro2010; Abderrahaman-Elena et al. Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019) and permeable substrates (Kuwata & Suga Reference Kuwata and Suga2017). However, this effect would likely be smeared out over canopies with irregularly aligned tips. The height over which the fluctuations decay within the canopy and the magnitude of the fluctuations in the core of the canopy appear to correlate with the element spacing. Note that this plateau in 
$u^{\prime }$ and 
$w^{\prime }$ within the canopy is asymptotic and requires a sufficiently large canopy depth to occur. Thus, this plateau is essentially absent for the canopy of S48, because of its low canopy height-to-spacing ratio, 
$h/s\approx 2$. The wall-normal fluctuations within the canopy do not exhibit this plateau and decay gradually below the canopy-tip plane to meet the impermeability condition at the canopy base. Let us also note here that the element-induced flow accounts for less than 30 % of the magnitude of the streamwise velocity fluctuations and less than 10 % of the cross-velocity fluctuations within the canopy, which is consistent with the observations of Poggi & Katul (Reference Poggi and Katul2008). This implies that the velocity fluctuations deep within the canopy result mainly from the penetration of the overlying, element-incoherent velocity fluctuations. This will be discussed further in § 3.3.
Although, as discussed in the preceding paragraphs, the element spacing has a leading-order effect on the fluctuating flow, their height, 
$h$, also plays a secondary role. In order to assess the effect of height, we consider the canopies of family H, which have fixed element width and spacing, but different canopy heights. As noted previously, the differences between the element-induced fluctuations for the fixed-spacing canopies of family H are negligible. However, we do observe changes in the full r.m.s. velocity fluctuations for these cases, implying that the height affects the element-incoherent flow. Above the canopy tips, we observe a decrease in 
$u^{\prime }$ and an increase in 
$v^{\prime }$ and 
$w^{\prime }$ with increasing canopy height, similar to the effect of increasing element spacing, as shown in figure 7(b,e,h). Within the canopy, 
$u^{\prime }$ and 
$w^{\prime }$ for all the cases collapse to the same curves, only departing to meet the no-slip condition at the canopy base. The corresponding magnitude of 
$v^{\prime }$ within the canopy, however, increases with canopy height up to 
$h/s\approx 6$, and saturates for 
$h/s\gtrsim 6$. This saturation is also observed for the effect of the canopy on the flow in general as illustrated in figures 7(b,e,h) and 9(b,e), which show that the velocity fluctuations, Reynolds shear stresses and the mean velocity profiles for cases H96 and H128, with 
$h/s\approx 6{-}8$ are essentially the same. The changes in the element-incoherent flow observed for different element heights likely result from a modulation of the Kelvin–Helmholtz-like instability, discussed in § 3.3, which is essentially independent of the element-induced flow.
Premultiplied spectral energy densities of the wall-normal velocity, 
$k_{x}k_{z}E_{vv}$, at height 
$y^{+}\approx 15$, normalised by their respective r.m.s. values. The line contours represent (a–e) cases S10 to S48; (f–j) cases H16 to H128; and (k–o) cases G10 to G100. The shaded contours represent the smooth-wall case, SC. The contours are in increments of 0.06 for all the cases. The vertical lines mark the most amplified wavelength predicted by linear stability analysis, discussed in § 4; ——, DNS mean profiles without drag on fluctuations; 
$\cdots \cdots$, DNS mean profiles with drag on fluctuations; - - -, synthesised mean profiles without drag on fluctuations.
3.3 Effect of canopy parameters on the shear-layer instability
The variations observed in the velocity fluctuations for the fixed-element-spacing simulations, discussed above, may result from the growth of the Kelvin–Helmholtz-like, shear-layer instability, typically reported in dense canopy flows (Finnigan Reference Finnigan2000; Nepf Reference Nepf2012). In order to assess the presence of this instability in the flow, we compare the premultiplied spectral energy densities of the wall-normal velocity at 
$y^{+}\approx 15$ in figure 10(f–j). For case H16 we observe that the spectral energy densities of the fluctuations above the canopy are similar to those above smooth walls. As the height of the canopy is increased, we observe a progressive increase in the energy in long spanwise wavelengths, 
$\unicode[STIX]{x1D706}_{z}^{+}>100$, for a narrow range of streamwise wavelengths, 
$\unicode[STIX]{x1D706}_{x}^{+}\approx 150{-}250$. This range of streamwise wavelengths remains roughly constant for increasing canopy heights. Such a signature in the spectral energy densities has been previously associated with the presence of spanwise-coherent, Kelvin–Helmholtz-like instabilities over riblets (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011), transitional roughness (Abderrahaman-Elena et al. Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019) and permeable substrates (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). This signature in the spectral energy densities is also reflected in the instantaneous realisations of the wall-normal velocity, portrayed in figure 5, which show increased spanwise coherence with increasing element height. The shear-layer instability is known to generate strong wall-normal fluctuations and, hence, its signature is most clear in the wall-normal spectra (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). Ghisalberti (Reference Ghisalberti2009) and Nepf (Reference Nepf2012) concluded that canopies only exhibit a shear-layer instability if their height is larger than the wall-normal extent of the instability as otherwise the rollers would be constrained by the lack of canopy depth. Above a short canopy, like that of H16, the proximity of the impermeability condition at the base of the canopy would inhibit the instability by blocking the wall-normal fluctuations. Similarly, Huerre (Reference Huerre1983) and Healey (Reference Healey2009) showed that the confinement of a free-shear layer also results in stabilisation of the associated Kelvin–Helmholtz instability. Increasing the canopy height weakens this effect, leading to a stronger signature of the instability, observed in figure 10(f–j). This enhanced signature of the instability is likely responsible for the increase in the cross-velocity fluctuations for the canopies family H with increasing height, discussed in § 3.2. For 
$h/s>6$, the instability no longer perceives the canopy base and the effect of the height on the instability and, consequently, the velocity fluctuations, saturates. The Kelvin–Helmholtz-like instability has also been reported to cause an increase in the Reynolds shear stresses, with an associated increase in the friction drag, over surfaces such as riblets and permeable substrates (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). The increase and saturation of the Reynolds shear stresses with increasing canopy height can be observed in figure 9(e) for the canopies of family H, and is concurrent with the effect of the element height on the instability, discussed above. This trend in the Reynolds shear stress has a corresponding effect on the drag exerted on the overlying flow, which is illustrated by the downward shift in the mean velocity profiles, portrayed in figure 9(b). The above discussion suggests that the secondary effect that the height has on the full velocity fluctuations within and above the canopy is mainly through its influence on the Kelvin–Helmholtz-like instability.
The increase in intensity of the Kelvin–Helmholtz-like rollers with increasing canopy height also contributes to the increase in the wall-normal velocity fluctuations within the canopy observed in figure 7(e). This is demonstrated by the wall-normal spectral energy densities of the flow within the canopies, portrayed at 
$y^{+}\approx -10$ for all the canopies of family H, in figure 11(a–e). Note that in calculating the spectra for a region with solid obstacles, we have implicitly assumed that the obstacles are fluid regions with zero flow velocity. As discussed in the previous paragraph, for case H16 the instability is inhibited by the proximity of the canopy-base wall, and the flow above shows similarities to a smooth-wall flow. The energy density within the canopy at 
$y^{+}\approx -10$ for this case also shows some regions overlapping with the smooth-wall spectra, with additional energy in the wavelengths associated with the Kelvin–Helmholtz-like instability. The smooth-wall spectra displayed for reference are at 
$y^{+}\approx 1$, which is as low as possible while yielding a non-negligible energy, since no direct comparison with 
$y^{+}\approx -10$ is possible. We also observe some energy in the spanwise wavelength corresponding to the canopy spacing and a broad range of streamwise wavelengths. These regions can be attributed to the modulation of the element-induced flow by the larger scale fluctuations induced by the instability or the overlying turbulence (Abderrahaman-Elena et al. Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019). This suggests that the fluctuations within a short canopy result mainly from the penetration of the overlying turbulence, with additional contributions from the Kelvin–Helmholtz-like instability and the element-induced flow. As the canopy height is increased, and the instability becomes stronger, the deviations in the spectral energy densities from smooth-wall flow become more prominent. The fluctuations within the canopy in cases H32 to H128 arise mainly from large spanwise wavelengths, likely originating from the Kelvin–Helmholtz-like instability near the canopy-tip plane, along with a contribution of the modulated element-induced flow discussed above. The increasing signature of the instability within the canopy with increasing element height can also be observed in the instantaneous realisations of the wall-normal velocity at 
$y^{+}\approx -10$ portrayed in figure 12. The presence of large spanwise wavelengths deep within the canopy can also be noted for the canopies of family S, whose spectral energy densities and realisations of wall-normal velocity at 
$y^{+}\approx -40$ are portrayed in figures 11(f–j) and 12, respectively. This suggests that, in the present dense canopies, the background turbulence is not able to penetrate far below the canopy tips, and that the velocity fluctuations deep within originate mainly from the footprint of the Kelvin–Helmholtz-like rollers above.
Premultiplied spectral energy densities of the wall-normal velocity, 
$k_{x}k_{z}E_{vv}$, for (a–e) cases H16 to H128 at a height of 
$y^{+}\approx -10$; and (f–j) cases S10 to S48 at a height of 
$y^{+}\approx -40$. The contours are normalised by the r.m.s. values of their respective cases. The shaded contours are of the smooth-wall case, SC at a height of 
$y^{+}\approx 1$, for reference. The contours are in increments of 0.075 for all the cases.
Instantaneous realisations of the wall-normal velocity at 
$y^{+}=-10$ (a,c,e,g,i) and 
$y^{+}=-40$ (b,d,f,h,j), normalised by 
$u_{\unicode[STIX]{x1D70F}}$. From top to bottom, (a,c,e,g,i) represent cases H16 to H128; and (b,d,f,h,j), cases S10 to S48. From top to bottom, the clearest and darkest colours indicate intensities of 
$\pm (0.1,0.2,0.3,0.3,0.3)$ in (a,c,e,g,i) and 
$\pm (0.05,0.2,0.4,0.4,0.5)$ in (b,d,f,h,j).
It is also worth noting that even in canopies with small element spacings, such as that of case S10, although the fluctuations of the wall-parallel velocities decay rapidly below the canopy-tip plane, the fluctuations of the wall-normal velocity decay more slowly, as shown in figure 7. This is also the case for the velocity fluctuations of the canopies of family H, for which the wall-normal velocity fluctuations within the canopy require larger canopy heights to saturate compared to the wall-parallel fluctuations. The presence of wall-normal fluctuations deep within the canopy are a reflection of the canopy layout being able to obstruct the wall-normal flow less efficiently than the tangential flow. It will be shown in § 4 that, for the present canopies, the effective drag coefficient in the tangential directions can be up to three times larger than in the wall-normal direction. In the core region of a tall canopy, the only mechanism to inhibit the velocity fluctuations is the canopy drag. As the canopy geometries studied here exert more drag on the wall-parallel flow than the wall-normal flow, 
$u^{\prime }$ and 
$w^{\prime }$ decay faster than 
$v^{\prime }$ within the canopy.
The spacing between the canopy elements affects the excitation of the Kelvin–Helmholtz-like instability through its influence on the canopy drag. White & Nepf (Reference White and Nepf2007) and Sharma, Gómez-de-Segura & García-Mayoral (Reference Sharma, Gómez-de-Segura and García-Mayoral2017) have shown that the canopy drag governs the instability through two competing effects, the shear at the canopy tips and the canopy drag. A small element spacing results in a large drag within the canopy, which in turn results in a larger shear at the canopy tips that enhances the instability, but at the same time it also inhibits the fluctuations more strongly, which weakens the instability. To study this effect, we now compare the premultiplied spectral energy densities of the wall-normal velocities for the canopies of family S, which have a constant height and different element spacing. For the canopy with the smallest spacing, S10, the signature of the Kelvin–Helmholtz-like instability in the energy densities is weak and the distribution of energy in different wavelengths is similar to that over smooth walls, as shown in figure 10(a). This suggests that the large drag exerted on the fluctuations by this canopy inhibits the formation of the shear-layer instability. As the element spacing is increased, the drag on the fluctuations reduces, and there is a progressive increase in the energy in wavelengths associated with the Kelvin–Helmholtz-like instability for the canopies of S16–S48, as can be observed in figure 10(b–e). The change in element spacing also has an effect on the instability wavelength, which will be discussed later in this section. In addition to the increase in the energy associated with the instabilities, the increase in element spacing also results in a progressive decrease in the overlapping regions in the energy densities of the canopy and smooth-wall flows, with a reduction in the energy in wavelengths 
$\unicode[STIX]{x1D706}_{x}^{+}\gtrsim 700$, 
$\unicode[STIX]{x1D706}_{z}^{+}\approx 50{-}100$. If the element spacing was increased further, the Kelvin–Helmholtz-like instability would eventually weaken, and for sparse enough canopies the flow within would begin to resemble smooth-wall flow perturbed by the element-induced flow of the isolated canopy elements. Such sparse canopies are beyond the scope of the present work, but have been discussed in Sharma & García-Mayoral (Reference Sharma and García-Mayoral2020), as well as in the previous studies of Poggi et al. (Reference Poggi, Porporato, Ridolfi, Albertson and Katul2004), Pietri et al. (Reference Pietri, Petroff, Amielh and Anselmet2009) and Huang et al. (Reference Huang, Cassiani and Albertson2009).
Let us now focus on the self-similar canopy geometries of family G. Although these canopies have the same 
$\unicode[STIX]{x1D706}_{f}$, increasing the size of the canopy elements produces similar effects on the premultiplied spectral energy densities as the canopies of family S, discussed in the previous paragraph. For the densest canopy, G10, the spectral energy densities at 
$y^{+}\approx 15$ collapse with those over a smooth wall, as shown in figure 13(a–d). As the size of the canopy is increased, we observe a stronger signature of the Kelvin–Helmholtz-like instability in the energy densities, portrayed in figure 10(k–o). The associated increase in spanwise coherence in the flow can also be observed in the instantaneous realisations of the wall-normal velocity shown in figure 5. Note that here, we observe the combined effects of increasing canopy height, as in family H, and increasing spacing, as in family S, on the instability. So far we have mainly discussed the spectral energy densities of the wall-normal velocity fluctuations, as they have the strongest signature of the Kelvin–Helmholtz-like instability. For completeness, we now use the canopies of family G to illustrate the effect of increasing the canopy size on the spectral energy densities of the streamwise and spanwise fluctuations, and the Reynolds shear stresses, portrayed in figure 13. The distinct region in the wall-normal spectral energy densities associated with the Kelvin–Helmholtz-like instability is not so apparent in the energy densities of the other velocity fluctuations and the Reynolds shear stresses. Nevertheless, as the canopy size increases, we observe an increase in the energy in streamwise wavelengths associated with the instability, along with a general increase in the energy in shorter and wider wavelengths compared to smooth-wall flows. This is consistent with the gradual shortening and widening of the eddies observed in figure 5.
Premultiplied spectral energy densities at 
$y^{+}\approx 15$ normalised by their respective r.m.s. value. The line contours represent (a–d) case G10; (e–h) case G40; (i–l) case G100. The filled contours represent the smooth-wall case, SC. The contours in (a,e,i), (b,f,j), (c,g,k) and (d,h,l) are in increments of 0.075, 0.06, 0.07 and 0.1, respectively.
The results discussed in this section suggest that the growth of the Kelvin–Helmholtz-like instability depends on both the canopy height and the element spacing. The streamwise wavelength of the instability, however, seems to depend mainly on the element spacing. We observe that, in the canopies of family H the streamwise wavelength of the instability is roughly constant regardless of the canopy height, 
$\unicode[STIX]{x1D706}_{x}^{+}\approx 150$, as can be observed in figure 10(f–j). For the canopies of family S, however, the increase in the element spacing results in an increase in the streamwise wavelength of the instability from 
$\unicode[STIX]{x1D706}_{x}^{+}\approx 140$ for case S10 to 
$\unicode[STIX]{x1D706}_{x}^{+}\approx 280$ for case S48, as shown in figure 10(a–e). Similarly, for the canopies of family G there is a progressive increase in the streamwise wavelengths associated with the instability with the increase in canopy size. The streamwise wavelength of the Kelvin–Helmholtz-like instability is determined by the shear length, typically defined in the literature as 
$L_{s}=U/(\text{d}U/\text{d}y)$ calculated at the canopy-tip plane (Raupach et al. Reference Raupach, Finnigan and Brunet1996; Finnigan Reference Finnigan2000; Nepf Reference Nepf2012). Previous studies have shown that the shear length 
$L_{s}$ in canopy flows can be determined by the effective streamwise canopy drag coefficient (Finnigan Reference Finnigan2000; Nepf et al. Reference Nepf, Ghisalberti, White and Murphy2007; Ghisalberti Reference Ghisalberti2009; Nepf Reference Nepf2012). Intuitively, in tall, dense canopies, we would expect this drag coefficient, and by extension the shear-layer thickness, to be a function of the element spacing. A dependence of the canopy shear-layer thickness on the element spacing was also observed by Novak et al. (Reference Novak, Warland, Orchansky, Ketler and Green2000) in their study of natural canopy flows. Therefore, the canopies of family H, which have a constant element spacing, have similar mean drag coefficients and, consequently, similar instability wavelengths, as observed in the spectral energy densities of the fixed-spacing canopies. For the canopies of families S and G, increasing the element spacing decreases the canopy drag coefficient, thereby resulting in the larger wavelengths observed for the shear-layer eddies. The effect of the canopy spacing on the drag coefficients and the instability wavelengths will be discussed further in § 4.
Finally, let us discuss the effect of the Reynolds number on the Kelvin–Helmholtz-like instability. It was shown in § 2.1, that the turbulent fluctuations over dense canopies scale in friction units, and therefore similar results are obtained when simulating canopies with the same height and spacing in friction units at different Reynolds numbers. It can be observed in figure 4(b) that the signature of this instability in simulations 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$ is essentially the same, and that the associated streamwise wavelength for both cases is roughly 
$\unicode[STIX]{x1D706}_{x}^{+}\approx 150$. As discussed above, the wavelength and amplification of the instability are governed by the shear at the canopy tips. As the canopy parameters for cases 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$ are kept constant in friction units, we can also expect the shear at the canopy tips to also be similar. Therefore, the instability characteristics for both these canopies are essentially the same when scaled in friction units. We have observed similar behaviours for Kelvin–Helmholtz-like instabilities originating over riblets (García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2012) and permeable substrates (Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019).
4 Linear analysis of Kelvin–Helmholtz-like instabilities
The results from DNS discussed in § 3 show that the flow in the region near the canopy-tip plane can be dominated by the presence of spanwise-coherent structures originating from a Kelvin–Helmholtz-like instability. This instability can be captured by a two-dimensional, mean-flow linear stability analysis, even in turbulent flows (Jimenez et al. Reference Jimenez, Uhlmann, Pinelli and Kawahara2001; White & Nepf Reference White and Nepf2007; García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011; Zampogna et al. Reference Zampogna, Pluvinage, Kourta and Bottaro2016; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). In this section, we discuss the methodology and results from such an analysis conducted on the velocity profiles obtained from the DNS. As the Kelvin–Helmholtz-like instability is an inviscid phenomenon, several of the studies just cited use an inviscid analysis to capture it. The inclusion of viscosity, however, inhibits the growth of smaller wavelengths in the flow, and consequently, results in the most amplified wavelength being slightly larger compared to that of an inviscid analysis (Jimenez et al. Reference Jimenez, Uhlmann, Pinelli and Kawahara2001; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). In this section, we present the results only from viscous analysis. The results from an inviscid analysis are presented in appendix B for reference. In addition, we show that some of the key features of this instability can be captured by linear analysis performed on velocity profiles modelled a priori, which would not require any information from the DNS.
For the purpose of the stability analysis, we model the effect of the canopy using a drag force in the Navier–Stokes equations, which results in the following governing equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D708}C_{i}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$ where 
$C_{i}$ is the effective canopy drag coefficient in each 
$i$th direction, has dimensions of inverse length squared – being essentially the inverse of a permeability – and is assumed to be homogeneous over the entire canopy region, as in Singh et al. (Reference Singh, Bandi, Mahadevan and Mandre2016). Given the density of the canopies considered, with maximum spacings 
$s^{+}=O(10)$, we assume that inertial effects in the flow deep within the canopy are small and can be neglected. In addition, the element width of the canopies is also small, 
$w^{+}\approx 1{-}25$. For such canopies, the canopy drag can be assumed to depend linearly on the velocity (Tanino, Nepf & Kulis Reference Tanino, Nepf and Kulis2005; Tanino & Nepf Reference Tanino and Nepf2008).
In the core of the canopy, away from the shear effects at the canopy base and top, the mean momentum equation would reduce to a balance between the canopy drag and the mean pressure gradient
$$\begin{eqnarray}\unicode[STIX]{x1D708}C_{x}U=-\frac{\text{d}P}{\text{d}x}.\end{eqnarray}$$ Equation (4.3) is essentially Darcy’s equation for flow within permeable substrates (Darcy Reference Darcy1856), and has been used by Zampogna & Bottaro (Reference Zampogna and Bottaro2016) to model flow deep within densely packed, rigid fibres. The streamwise drag coefficient, 
$C_{x}$, can be obtained by substituting the values of 
$U$ and 
$\text{d}P/\text{d}x$ obtained from the DNS into (4.3). From dimensional arguments, equation (4.3) predicts that the drag coefficient would scale as 
$C_{x}\sim 1/s^{2}$. This scaling is demonstrated in figure 14(a), which suggests that (4.3) provides a reasonable approximation for the flow deep within the present canopies, excluding the sparsest canopy S48. Although we can expect the flow within the canopy to be Darcy-like in the wall-normal direction as well, we cannot use the DNS results to obtain 
$C_{y}$, as there is no mean flow in this direction. In order to obtain 
$C_{y}$, we consider separately the Stokes flow along infinitely long canopy elements driven by a constant pressure gradient. The equation for such flow is 
$\unicode[STIX]{x1D708}(\unicode[STIX]{x2202}_{x}^{2}+\unicode[STIX]{x2202}_{z}^{2})v=\text{d}P/\text{d}y$. The wall-normal drag coefficient is then obtained as
$$\begin{eqnarray}\unicode[STIX]{x1D708}C_{y}\langle v\rangle =-\frac{\text{d}P}{\text{d}y},\end{eqnarray}$$ where the angled brackets represent a spatial average. The estimated values of 
$C_{y}$ are portrayed in figure 14(b) for reference. It may be noted that the ratio of the streamwise to wall-normal drag coefficients for the present canopies is 
$C_{x}/C_{y}\approx 2{-}3$, which shows that the streamwise flow is more obstructed than the wall-normal flow for the layouts considered. It is worth noting here that the canopy drag coefficients, 
$C_{x}$ and 
$C_{y}$, and the ratio between them also depends significantly on the plan view arrangement and the resulting porosity of the canopy (Van der Westhuizen & Du Plessis Reference Van der Westhuizen and Du Plessis1996; Zampogna & Bottaro Reference Zampogna and Bottaro2016). This is evidenced by the different ratios of the drag length scale and the element spacing for the canopies of families S and G, portrayed in figure 14, which have width-to-spacing ratios, 
$w/s=1/2$ and 
$2/9$, respectively. The dependence of the 
$C_{x}$ on 
$w/s$ can also be predicted using two-dimensional Stokes flow simulations, as shown in figure 14(a). For more complicated canopy arrangements, such as staggered or random, we would expect the drag to be a function of the planar layout of the elements.
Variation of the length scales derived from the (a) streamwise and (b) wall-normal canopy drag coefficients for different element spacings. The symbols represent, ▫, cases of S; 
$+$, cases of H; and ○, cases of G. The colours from red to blue represent cases S10 to S48, H16 to H128 and G10 to G100. The symbols in (a) are values obtained from the DNS, and the dashed and solid lines are predictions from two-dimensional Stokes-flow simulations. Both the symbols and the lines in (b) are obtained from Stokes-flow simulations.
In order to conduct the stability analysis, we linearise the equations (4.1) and (4.2) around the mean flow, 
$U(y)$, yielding
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+U\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+vU^{\prime }=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}u-\unicode[STIX]{x1D708}C_{x}u, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}t}+U\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}x}=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}v-\unicode[STIX]{x1D708}C_{y}v, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}y}=0. & \displaystyle\end{eqnarray}$$These equations are used to obtain a modified Orr–Sommerfeld equation (Drazin & Reid Reference Drazin and Reid2004; White & Nepf Reference White and Nepf2007; Singh et al. Reference Singh, Bandi, Mahadevan and Mandre2016; Zampogna et al. Reference Zampogna, Pluvinage, Kourta and Bottaro2016),
$$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+U\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D708}C_{y}\right)\unicode[STIX]{x1D6FB}^{2}v-\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{4}v=U^{\prime \prime }\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D708}(C_{x}-C_{y})\frac{\unicode[STIX]{x2202}^{2}v}{\unicode[STIX]{x2202}y^{2}}. & & \displaystyle\end{eqnarray}$$ Assuming wave-like solutions of the form 
$v=\widetilde{v}\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x-\unicode[STIX]{x1D714}t)}$, equation (4.8) reduces to the eigenvalue problem
$$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D6FC}U-\text{i}\unicode[STIX]{x1D708}C_{y})(\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2})\widetilde{v}-\unicode[STIX]{x1D6FC}U^{\prime \prime }\widetilde{v}-\text{i}\unicode[STIX]{x1D708}(C_{x}-C_{y})D^{2}\widetilde{v}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\text{i}\unicode[STIX]{x1D708}(\text{D}^{4}-2\unicode[STIX]{x1D6FC}^{2}\text{D}^{2}+\unicode[STIX]{x1D6FC}^{4})\widetilde{v}=\unicode[STIX]{x1D714}(\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2})\widetilde{v},\end{eqnarray}$$ where the prime superscript denotes differentiation with respect to 
$y$, and 
$\text{D}$ represents the operator 
$\text{d}/\text{d}y$. Equation (4.9) is solved to obtain the complex frequency, 
$\unicode[STIX]{x1D714}$, for real values of the streamwise wavenumber, 
$\unicode[STIX]{x1D6FC}$, subject to no-slip and impermeability boundary conditions at the top and bottom walls. The instability is then amplified for positive values of the imaginary part of 
$\unicode[STIX]{x1D714}$.
Growth rates of different perturbation wavelengths obtained from the stability analysis performed on (a–c) mean profiles obtained from the DNS, with drag on the perturbations included in the stability analysis; (d–f) mean profiles obtained from DNS, with no drag on the perturbations; and (g–i) mean velocity profiles obtained using (4.10), with no drag on the perturbations. The lines from red to blue, indicated by the direction of the arrows, represent (a,d,g) cases S10 to S48; (b,e,h) cases H16 to H128; and (c,f,i) cases G10 to G100.
Contours of the stream function for the most amplified mode for case H96 obtained from the stability analysis (a) with drag and (b) without drag on the perturbations. The blue and red lines correspond to clockwise and counter-clockwise rotation, respectively.
Most amplified instability wavelengths observed in the DNS and predicted by the stability analysis, scaled in friction units. The column labelled ‘DNS’ lists the approximate streamwise wavelength associated with the instability in the wall-normal spectra portrayed in figure 10. Here, 
$\text{SA}_{C0}$, most amplified wavelengths from stability analysis on DNS mean profiles without drag on fluctuations; 
$\text{SAM}_{C0}$, on synthesised velocity profiles without drag on fluctuations; and SA, on DNS mean profiles with drag on fluctuations.
The growth rates for different perturbation wavelengths are portrayed in figure 15(a–c), and the wavelengths with the highest growth rates are summarised in table 2. The most amplified wavelengths predicted by the stability analysis only match those observed in the DNS for canopies with high values of 
$\unicode[STIX]{x1D6FF}/h$. The wavelengths predicted for cases H16, H32, G10, G20 and G40 show reasonable agreement with those observed in the DNS. For canopies with larger heights, however, the analysis predicts wavelengths larger than those observed in the DNS. For the fixed-spacing canopies of family H, the predicted instability wavelength also increases with increasing canopy height, whereas the DNS show that the instability wavelength for these cases is essentially independent of the height. The contours of the instability stream function for case H96 for the most amplified wavelength, 
$\unicode[STIX]{x1D706}_{x}^{+}\approx 385$, portrayed in figure 16(a), show that it has a large wall-normal span, extending up to 
$y^{+}\approx 120$. Such an instability was also reported by Singh et al. (Reference Singh, Bandi, Mahadevan and Mandre2016), who performed stability analyses similar to the one conducted here, except that the canopy was represented by a drag force depending quadratically on the velocity. Singh et al. (Reference Singh, Bandi, Mahadevan and Mandre2016) noted that their analysis predicted two instability modes, one similar to the Kelvin–Helmholtz instability and another originating from the canopy drag included in the analysis. They only considered canopies with low 
$\unicode[STIX]{x1D6FF}/h$, and observed that the second instability mode, similar to the large-wavelength modes obtained from the present stability analysis, was dominant for canopies with high drag and spanned the entire height of the channel.
It is worth noting here, that in the region near the interface between the canopy and the free-flow, the assumption of a constant drag coefficient given by (4.3) would no longer be valid, as shear and advective effects become stronger. We have conducted some exploratory analyses accounting for this variation in the drag coefficient and these do not provide improved estimates for the instability wavelength compared to the results of the constant-drag analysis presented here. In the present analysis, we have also assumed that the drag coefficient experienced by the perturbations is the same as that experienced by the mean flow. However, we have recently reported for sparser canopies that different wavelengths in the flow can perceive drag coefficients different from that for the mean flow (Sharma & García-Mayoral Reference Sharma and García-Mayoral2020). In such a case, the drag coefficient would have to be calculated on a mode-by-mode basis for the different wavelengths. A wavelength-dependent drag coefficient would lead to the drag for a given wavelength being dependent on a convolution from all other wavelengths. Such an analysis, however, is beyond the scope of the present work. An alternative approach would be to model the canopy as a permeable substrate, which naturally yields wavelength-dependent equations for the flow within (Zampogna et al. Reference Zampogna, Pluvinage, Kourta and Bottaro2016; Abderrahaman-Elena & García-Mayoral Reference Abderrahaman-Elena and García-Mayoral2017; Sharma et al. Reference Sharma, Gómez-de-Segura and García-Mayoral2017; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). In order to illustrate that applying a constant drag coefficient on the perturbations may be a coarse assumption, we also present results from stability analyses with no drag on the perturbations. This is also a rather coarse assumption, but we observe that excluding the drag on the perturbations in the stability analysis yields better estimates for the instability wavelengths observed in the DNS, as shown in figure 10. For the canopies of family H, the stability analysis without drag on the fluctuations shows that the most amplified wavelength does not vary significantly with the canopy height. For the canopies of families S and G, this analysis shows an increase in the most amplified wavelength with increasing element spacing, owing to the increase in the shear-layer thickness. The results from this analysis are portrayed in figure 15(d–f), and the most amplified wavelength for each case is listed in table 2. While neglecting the drag acting on the fluctuations yields better estimates for the most amplified wavelengths for canopies with small spacings, the predictions for larger spacings differ by up to a factor of two from the direct numerical simulation observations. This is likely due to the assumption that the mean flow is homogeneous in the tangential directions, implicit in the stability analysis, which breaks down for such cases. We have not observed any significant signature of the Kelvin–Helmholtz-like instability for the sparser canopies studied in Sharma & García-Mayoral (Reference Sharma and García-Mayoral2020) despite the presence of an inflection in the mean velocity profiles. There may also be some distortion of the instability by the ambient turbulent fluctuations in the DNS (Rogers & Moser Reference Rogers and Moser1994; Raupach et al. Reference Raupach, Finnigan and Brunet1996).
Growth rate for different perturbation wavelengths from the stability analysis for cases ——, 
$\text{H}32_{180}$; and - - -, 
$\text{H}32_{400}$; (a) with drag and (b) without drag on the perturbations.
We have also performed stability analyses on the cases with different Reynolds numbers, 
$\text{H}32_{180}$ and 
$\text{H}32_{400}$. These analyses predict similar instability wavelengths and growth rates in viscous units, as shown in figure 17, which is consistent with the observations in the DNS, discussed in § 3.3. These results emphasize that the strength of inflection in the mean-velocity profile and the shear-layer thickness for both these cases is similar when scaled in friction units and, hence, so is their effect on the instability.
The results obtained from the DNS and the stability analysis suggest that there is a dependence on the element spacing of the most amplified wavelength, related to the effect of the spacing on the shear-layer thickness. The usual definition of the shear-layer thickness, 
$L_{s}=U/(\text{d}U/\text{d}y)$, misses the contribution of the part of the shear-layer above the canopy. Regarding the latter, García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011) studied the formation of Kelvin–Helmholtz-like instabilities over riblets, and noted that the shear-layer thickness above was given by the height at which the vorticity gradient, 
$\text{d}^{2}U/\text{d}y^{2}$, concentrated. In smooth-wall flows, this height is roughly 
$y_{c}^{+}\approx 5{-}10$. For the present cases, we observe that the instability wavelengths predicted by the stability analysis correlate well with the full shear length 
$L_{s}+y_{c}$, if we take 
$y_{c}^{+}\approx 5$, as shown in figure 18(a). This suggests that the shear-layer semi-thickness above the canopies, 
$y_{c}$, is roughly constant for most of the geometries considered here, and remains close to the smooth-wall value, while the semi-thickness below has the standard form 
$L_{s}=U/(\text{d}U/\text{d}y)$, measured at the canopy tips plane. The only notable deviation is for the sparsest canopy studied, S48. For canopies with large element spacings, we observe that the peak in 
$\text{d}^{2}U/\text{d}y^{2}$ moves closer to the canopy-tip plane, so 
$y_{c}^{+}=5$ may no longer be a reasonable approximation for the shear-layer semi-thickness above. Regarding the height of the shear layer within the canopy, 
$L_{s}$, we observe that it is set by the mean canopy drag coefficient, 
$L_{s}\propto \sqrt{1/C_{x}}$, as also noted in the studies of aquatic canopy flows by Nepf et al. (Reference Nepf, Ghisalberti, White and Murphy2007) and White & Nepf (Reference White and Nepf2007). The drag coefficient on the mean flow, in turn, depends on the element spacing and the width-to-spacing ratio, as shown in figures 18(b) and (c). The correlation of 
$L_{s}$ with 
$s$, therefore, explains the dependence of the most amplified wavelength on the element spacing observed in the DNS and the stability analysis.
(a) Instability wavelength, 
$\unicode[STIX]{x1D706}_{x}^{+}$, obtained from the linear stability analysis versus the total shear length, 
$L_{s}^{+}+y_{c}^{+}$; (b) shear length, 
$L_{s}^{+}$, versus the drag length scale; (c) shear length versus the element spacing. ▫, family S; 
$+$, family H; ○, family G. The colours from red to blue represent cases S10 to S48, H16 to H128 and G10 to G100. In (a) and (b), the solid lines are linear regressions with slopes 0.06 and 1.36, respectively. In (c), the solid and dashed lines are linear regressions with slopes 0.14 and 0.3, respectively.
4.1 Analysis on modelled velocity profiles
In this section we introduce a simple model for the mean velocity profile in dense canopy flows and discuss the results from their stability analysis. As we only consider canopies with small element spacings, 
$s^{+}=O(10)$, the magnitude of inertial effects within the canopies are also small and are thus neglected in the model. The results discussed in § 3 also suggest that, for very dense canopies, turbulence and, consequently, the Reynolds shear stresses do not penetrate within (Nepf et al. Reference Nepf, Ghisalberti, White and Murphy2007), and are smooth-wall-like above the canopy-tip plane. The mean velocity above the canopy could then be modelled using a smooth-wall eddy viscosity, with the canopy-tip plane acting as the location of the smooth-wall (Jimenez et al. Reference Jimenez, Uhlmann, Pinelli and Kawahara2001; García-Mayoral & Jiménez Reference García-Mayoral and Jiménez2011; Gómez-de-Segura, Sharma & García-Mayoral Reference Gómez-de-Segura, Sharma and García-Mayoral2018; Gómez-de-Segura & García-Mayoral Reference Gómez-de-Segura and García-Mayoral2019). The equation for the mean velocity can then be written as
$$\begin{eqnarray}\frac{\text{d}}{\text{d}y}\left([\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{T}(y)]\frac{\text{d}U}{\text{d}y}\right)-\unicode[STIX]{x1D708}C_{x}(y)U-\frac{\text{d}P}{\text{d}x}=0,\end{eqnarray}$$ where 
$C_{x}(y)$ is the average streamwise canopy drag coefficient, which is assumed constant within the canopy and zero outside, and 
$\unicode[STIX]{x1D708}_{T}(y)$ is the height-dependent eddy viscosity proposed by Cess (Reference Cess1958) to approximate turbulent smooth-channel flow, and is non-zero only outside the canopy. The drag coefficients, 
$C_{x}$, used to obtain the velocity profiles are those given by (4.3) and portrayed in figure 14(a). These have been obtained using the data from the DNS but can also be obtained from Stokes-flow simulations as shown in figure 14(a), which are significantly less computationally intensive.
The most amplified wavelengths predicted by the stability analysis conducted on these modelled velocity profiles, with no drag applied on the fluctuations, are in reasonable agreement with those obtained from the same no-drag analysis on profiles obtained from the DNS. The growth rates predicted are portrayed in figure 15(g–i). The wavelengths with maximum growth rates are also summarised in table 2. We have also conducted stability analyses on these modelled velocity profiles including the effect of the eddy viscosity. The results are portrayed in appendix B, and they are essentially the same as the ones obtained using molecular viscosity alone, apart from a slight reduction in the instability growth rates, which suggests that although 
$\unicode[STIX]{x1D708}_{T}$ is important for setting the shape of the mean velocity profile, its effect on the fluctuations is not significant. This is likely because the Kelvin–Helmholtz-like rollers occur near the canopy tips, where 
$\unicode[STIX]{x1D708}_{T}$ is small and the molecular viscosity, 
$\unicode[STIX]{x1D708}$, dominates. It is worth noting that even though this model is able to capture the instability wavelength, the velocity profiles obtained using this model do not match those from the DNS, apart from those of S10 and G10. This is most likely due to our assumption that the turbulent stresses do not penetrate within the canopy and remain smooth-like, which fails as the element spacing is increased. As discussed previously, the wavelength of the instability is set by the shear length. The shear-layer semi thickness within the canopy, 
$L_{s}$, is set by the canopy drag coefficient, 
$C_{x}$. As this drag coefficient is the same both from the DNS and for the modelled velocity profiles, we expect 
$L_{s}$ to be similar as well. The shear length above the canopy, however, could differ, as the profiles from DNS would include the effect of the turbulent stresses penetrating into the canopy and deviating from their smooth-wall values, while the modelled velocity profiles do not. The similarity in the instability wavelengths between these analyses therefore suggest that, for most of the dense canopies considered in this work, turbulence is essentially precluded from penetrating into the canopy, and that the shear length above does not vary significantly from its smooth-wall value.
5 Conclusions
In the present work, we have examined the effect of the canopy layout on turbulent flows over canopies of densely packed filaments of small size. Three families of simulations have been conducted, the first with the element height in friction units fixed, the second with the element spacing fixed, and the third with the height-to-spacing ratio fixed. The layouts considered had height-to-spacing ratios greater than one, and elements spacings in the range 
$s^{+}\approx 3{-}50$. The penetration of turbulent fluctuations within such canopies was limited by their small element spacings. Consequently, the height of the roughness sublayer was also determined by the element spacing, rather than their height, extending up to 
$y\approx 2{-}3s$ above the canopy tips. The canopy drag coefficient was also found to be determined by the element spacing, 
$C_{x}\sim 1/s^{2}$. Canopies with small spacings were, therefore, found to suppress the velocity fluctuations within them owing to the large drag exerted, and the fluctuations became more intense as the spacing increased. The intensity of the characteristic Kelvin–Helmholtz-like instability over canopies was observed to be governed by two competing effects resulting from the canopy drag, the inflection at the canopy tips and the drag on the fluctuations. Canopies with large drag had a large shear at the canopy tips, and thus a stronger inflection, which enhanced the instability, but also exerted a large drag on the velocity fluctuations, which suppressed the instability. The instability was found to be inhibited in canopies with 
$s^{+}\lesssim 10$ and, for the range of canopy spacings considered here, a stronger signature of the instability was observed as the spacing was increased. We also showed that the main contribution to the velocity fluctuations deep within the canopy was the footprint of the Kelvin–Helmholtz-like instability, and that the contribution of the element-induced dispersive flow was negligible. Short canopies with 
$h/s\sim 1$ were also found to inhibit the instability, owing to the blocking effect of the wall at the canopy base. For height to spacing ratios 
$h/s\gtrsim 6$, the instability was no longer influenced by the bottom wall, and the effect of the canopy height on the flow within and above the canopy saturated. Increasing the canopy height for a fixed spacing did not change the element-induced velocity fluctuations, and instead affected the surrounding flow through the influence of height on the instability.
Linear stability analysis conducted on the mean velocity profiles obtained from the DNS is able to capture the approximate wavelength of the instability observed in the DNS for canopies with small element spacings. The analysis fails for larger element spacings, for which the assumption of the flow perceiving the canopy in a homogenised fashion breaks down. We showed that the shear-layer thickness, which determines the instability wavelength, has two components, one within the canopy and the other above. The latter is set by the height above the canopy tips at which the vorticity gradient concentrates, and is essentially constant for the present canopies, 
$y_{c}^{+}\approx 5$. The shear-layer thickness within the canopy follows the conventional definition, 
$L_{s}=U/(\text{d}U/\text{d}y)$, and is determined by the canopy drag, thus depending linearly on the canopy spacing. We have also proposed a simplified model to capture the most amplified instability wavelength over dense canopies. The model assumes that the turbulence above the canopy does not penetrate within and remains smooth-wall-like, and uses the mean streamwise drag coefficient of the canopies to synthesise an approximate mean-flow profile. The stability analysis conducted using these synthesised profiles yields similar results to those conducted using the mean profiles from DNS.
Acknowledgements
A.S. was supported by an award from the Cambridge Commonwealth, European and International Trust. Computational resources were provided by the ‘Cambridge Service for Data Driven Discovery’ operated by the University of Cambridge Research Computing Service and funded by EPSRC Tier-2 grant EP/P020259/1.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Immersed-boundary method and validation
For the present work, a modified version of the immersed-boundary algorithm proposed by García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011) is used. The algorithm of García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011) was based on a direct-forcing approach, which applies a body force within the immersed-boundary points to drive the velocity at these points to zero (Mittal & Iaccarino Reference Mittal and Iaccarino2005). The condition to implement at the points within the canopy elements is
$$\begin{eqnarray}\frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\unicode[STIX]{x0394}t}=\frac{-\boldsymbol{u}^{n}}{\unicode[STIX]{x0394}t}.\end{eqnarray}$$Following García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011), this condition can be approximated by modifying the right-hand side of (2.1)
$$\begin{eqnarray}\left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}\right]\boldsymbol{u}_{k}^{n}=-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}\boldsymbol{u}_{k-1}^{n}.\end{eqnarray}$$Velocity profiles obtained using the two immersed-boundary algorithms described in appendix A, after one time step, starting from random initial conditions. Here, (a,b) show results obtained from the algorithm utilised by García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011) and Abderrahaman-Elena et al. (Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019), given by (A 2), and (c,d) those from the present algorithm, using (A 5). The shaded regions mark the location of the solid obstacles. The same data are plotted in the (a,c) and (b,d) columns, except that the right column portrays the velocities in a logarithmic scale.
The original code of García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011) used a collocated grid for the wall-normal coordinate and was extended by Fairhall & García-Mayoral (Reference Fairhall and García-Mayoral2018) to employ a staggered grid. Fairhall & García-Mayoral (Reference Fairhall and García-Mayoral2018) also split the Laplacian operator on the left-hand side of (2.1) into its wall-parallel and wall-normal components following Kim & Moin (Reference Kim and Moin1985)
$$\begin{eqnarray}\left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}\right]\boldsymbol{u}\approx \left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}_{xz}\right]\left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}_{y}\right]\boldsymbol{u},\end{eqnarray}$$ where 
$L_{xz}$ includes the wall-parallel components of 
$L$, and 
$L_{y}$ the wall-normal one. Splitting the Laplacian in this manner still retains the second-order temporal accuracy of the code (Kim & Moin Reference Kim and Moin1985). Equation (2.1) can then be written as
$$\begin{eqnarray}\left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}_{y}\right]\boldsymbol{u}=\left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}_{xz}\right]^{-1}\text{RHS},\end{eqnarray}$$where RHS refers to the right-hand side of (2.1). In the present work, we implement a modified version of the immersed-boundary algorithm used by García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011) into the above algorithm, which offers an improvement in the accuracy for the velocities within the immersed-boundary regions. This implementation is summarised below. The right-hand side of (A 4) is then transformed to physical space, and modified to satisfy the following conditions within the immersed-boundary points
$$\begin{eqnarray}\left[\unicode[STIX]{x1D644}-\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}_{y}\right]\boldsymbol{u}_{k}^{n}=\left\{\begin{array}{@{}ll@{}}\displaystyle -\unicode[STIX]{x0394}t\frac{\unicode[STIX]{x1D6FD}_{k}}{Re}\text{L}_{y}\boldsymbol{u}_{k-1}^{n}, & y=y_{tips},\\ 0, & y\neq y_{tips},\end{array}\right.\end{eqnarray}$$ where 
$y_{tips}$ denotes the wall-normal plane at the canopy tips. At the interfaces, the condition imposed by (A 5) yields 
$\boldsymbol{u}_{k}^{n}=O(\unicode[STIX]{x0394}t^{2})$, which is of the order of the temporal discretisation error of the code. As a staggered grid is used, the wall-normal grid points for the streamwise velocities are offset by half a grid spacing from those of the wall-normal velocity. The element tips are aligned with the grid for the streamwise velocity. For the wall-normal velocity, the interface condition is set at the grid point just below the canopy-tip plane, which enforces near-zero wall-normal velocity at the canopy tips through continuity. Away from the interfaces, within the immersed boundary region, the condition set by (A 5) results in an exponential decay of the velocity in the wall-normal direction from its 
$O(\unicode[STIX]{x0394}t^{2})$ value at the interface. To illustrate this decay, the results from a simple, one-dimensional implementation of this algorithm are shown in figure 19. The velocity is assumed to vary only in the wall-normal direction. For these test cases, random initial conditions are used to mimic turbulent fluctuations in the flow. The velocity fields obtained for these cases, after one time step, using the immersed boundary conditions of (A 5) are compared to those obtained using (A 2) in figure 19.
Although both algorithms result in small velocities within the solid regions, the implementation proposed here results in a smoother and much faster decay of the velocity within the solid. In the DNS code, however, the velocity correction step introduces an error of order 
$\unicode[STIX]{x0394}t^{2}$ in all the immersed-boundary points. Even so, in experience it was observed that the proposed algorithm is a more stable numerical implementation of the immersed boundaries compared to the one proposed by García-Mayoral & Jiménez (Reference García-Mayoral and Jiménez2011). This is likely due to the present method not generating sharp gradients in the velocity field within the solid obstacles at the pressure calculation step. The velocity within the canopy elements, or the permeability error, in the DNS is observed to be less than 
$0.1u_{\unicode[STIX]{x1D70F}}$, for all the conducted simulations, and is much smaller than velocity in the ‘fluid’ points surrounding the elements. This is illustrated in the instantaneous realisations of the velocity fields from one of the simulations are portrayed in figure 20. For completeness, the wall-normal grid distribution for case S10 is portrayed in figure 21. In order to validate the implementation of the immersed boundaries, we have replicated the DNS for the collocated roughness elements with height 
$h^{+}\approx 12$ of Abderrahaman-Elena et al. (Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019). The mean velocity profiles, r.m.s. fluctuations and the Reynolds shear stresses obtained from the present simulations show good agreement with the results of Abderrahaman-Elena et al. (Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019), as shown in figure 22. We have also conducted a grid-dependence analysis, for which the results are portrayed in figure 23. Velocity r.m.s. fluctuations for cases S48 and G100 are shown for different wall-parallel resolutions.
Instantaneous realisations of the (a) streamwise, (b) wall-normal and (c) spanwise velocities in a plane passing through the middle of the canopy elements for case S48, scaled with the friction velocity 
$u_{\unicode[STIX]{x1D70F}}$. The clearest and darkest contours represent intensities of 
$\pm 0.1$, respectively.
Wall-normal grid distribution for the simulation of case S10. (a) Variation of the wall-normal coordinate and (b) grid resolution with an equispaced auxiliary variable j. Dashed lines mark the location of the canopy tip plane.
Root mean square velocity fluctuations, mean velocity and Reynolds shear stress profiles. The solid lines represent the results obtained from the present code, and the 
$+$ symbol represent the data of case C12 from Abderrahaman-Elena et al. (Reference Abderrahaman-Elena, Fairhall and García-Mayoral2019).
Root mean square velocity fluctuations and Reynolds shear stress profiles. The panels (a,c,e,g) represent case G100 and those in (b,d,f,h) represent case S48. The solid lines in (a,c,e,g) represent the results from using 27 points per spacing; the symbols ○, 18 points; and 
$+$, 9 points per spacing. The solid lines in (b,d,f,h) represent the results from using 36 points per spacing; the symbols ○, 24 points; and 
$+$, 12 points per spacing.
Growth rates of different perturbation wavelengths obtained from the stability analysis performed on (a–c) mean profiles obtained from the DNS, with drag on the perturbations included in the stability analysis; (d–f) mean profiles obtained from DNS, with no drag on the perturbations; and (g–i) mean velocity profiles obtained using equation (4.10), with no drag on the perturbations. ——, viscous analysis including molecular viscosity alone; - - -, inviscid analysis; 
$\cdots \cdots$, viscous analysis including an eddy viscosity. The colours from red to blue represent (a,d,g) cases S10 to S48; (b,e,h) cases H16 to H128; and (c,f,i) cases G10 to G100.
Appendix B. Comparison of inviscid and viscous stability analysis
Here, we provide the governing equations used to perform a stability analysis with a turbulent viscosity varying in the wall-normal direction, 
$\unicode[STIX]{x1D708}_{T}(y)$. The modified Orr–Sommerfeld equation is then given by
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6FC}U-\text{i}C_{y})(\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2})\widetilde{v}-\unicode[STIX]{x1D6FC}U^{\prime \prime }\widetilde{v}-\text{i}(C_{x}-C_{y})D^{2}\widetilde{v}+2\text{i}\unicode[STIX]{x1D708}_{T}^{\prime }(\text{D}^{3}-\unicode[STIX]{x1D6FC}^{2}\text{D})\widetilde{v} & \displaystyle \nonumber\\ \displaystyle & \displaystyle +\,\text{i}\unicode[STIX]{x1D708}_{T}(\text{D}^{4}-2\unicode[STIX]{x1D6FC}^{2}\text{D}^{2}+\unicode[STIX]{x1D6FC}^{4})\widetilde{v}+\text{i}\unicode[STIX]{x1D708}_{T}^{\prime \prime }(\text{D}^{2}+\unicode[STIX]{x1D6FC}^{2})\widetilde{v}=\unicode[STIX]{x1D714}(\text{D}^{2}-\unicode[STIX]{x1D6FC}^{2})\widetilde{v}. & \displaystyle\end{eqnarray}$$A similar equation, excluding the canopy drag terms, has also been used by Reynolds & Hussain (Reference Reynolds and Hussain1972), Del Alamo & Jiménez (Reference Del Alamo and Jiménez2006), Pujals et al. (Reference Pujals, García-Villalba, Cossu and Depardon2009) and Gómez-de-Segura et al. (Reference Gómez-de-Segura, Sharma and García-Mayoral2018). In figure 24(a–f), we compare the results obtained from viscous and inviscid analysis conducted using the velocity profiles from the DNS. In figure 24(g–i), we show the results from inviscid and viscous stability analyses, using both molecular and turbulent viscosities, performed using the velocity profiles obtained from (4.10).
Open access
