Turbulent flows over dense filament canopies

Turbulent flows over dense canopies consisting of rigid filaments of small size are investigated using direct numerical simulation. The effect of the height and spacing of the canopy elements on the overlying turbulence is studied. In these dense canopies, the element spacing is the key parameter determining the penetration of turbulent fluctuations within the canopy and the height of the roughness sublayer. The flow within and above the canopy is found to become independent of the filament height at large element height-to-spacing ratios, $h/s \gtrsim 6$. We also study the effect of the canopy parameters on the Kelvin$-$Helmholtz-like instability typically associated with dense-canopy flows. We find that the instability is inhibited in canopies with small element height-to-spacing ratios, $h/s \sim 1$, due to the blocking effect of the wall at the base of the canopy, while a stronger signature of the instability can be observed as the height is increased. Canopies with very small element spacings, $s^+ \lesssim 10$, also inhibit the instability, owing to the large drag they exert on the fluctuations. The streamwise wavelength of the instability depends on the canopy shear-layer thickness, and we show that, for the present dense canopies, this thickness also scales with the element spacing. Some of the effects of the canopy parameters on the instability can be captured by a linear stability analysis.


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 (2000), Belcher et al. (2012) and Nepf (2012). In engineering applications, filament canopies can, for instance, be used to enhance heat transfer (Fazu & Schwerdtfeger 1989;Bejan & Morega 1993) and for energy harvesting (McGarry & Knight 2011;Elahi et al. 2018). Depending on the geometry and spacing of their elements, canopies can be classified as sparse, dense or transitional (Nepf 2012). Dense canopies typically have small element spacings compared to the lengthscales 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. 2004;Nepf 2012;. Transitional, or intermediate, canopies would lie between these two regimes. Flows over dense canopies are characterised by the formation of Kelvin-Helmholtz-like, mixing-layer instabilities near the top of the canopy (Finnigan 2000;Nepf 2012). These instabilities originate from the inflection point in the mean velocity profile at the canopy-tip plane (Raupach et al. 1996). Kelvin-Helmholtz instabilities manifest as spanwise coherent rollers whose streamwise scale is determined by the shear-layer thickness (Michalke 1972;Brown & Roshko 1974). Ghisalberti & Nepf (2004) 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. 1996;Singh et al. 2016;Luminari et al. 2016). 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 et al. 2009;Bailey & Stoll 2016). The importance of this instability decreases as the element spacing is increased, and sparse canopies do not exhibit a notable signature (Poggi et al. 2004;Pietri et al. 2009;Huang et al. 2009;Sharma & García-Mayoral 2019). It could be expected, however, that the effect of increasing the canopy height at a fixed element spacing on this instability and the surrounding canopy flow eventually saturated. Such an effect was noted by Sadique et al. (2017), who examined the meanvelocity profiles in turbulent flows over high-aspect-ratio prismatic posts. Sadique et al. (2017) found that the mean-velocity profile over such geometries 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. (2017) observed the height of this active region to be related to the element width. A similar observation was also made by MacDonald et al. (2018), who performed direct numerical simulations of flows over spanwise-aligned bars. They found that the gap between the bars was the relevant lengthscale for the overlying flow. 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 produced by them.
In the present work, we study turbulent flows over dense canopies. The canopies consist of rigid, prismatic filaments with small element spacings and large height-tospacing ratios. The spacings range from s + ≈ 2.6 to 48 for the densest and sparsest canopy geometries considered, respectively. We assess the effect of the element height and spacing on the turbulence within the canopies and on the overlying flow. The canopy geometries considered here are dense enough to trigger the Kelvin-Helmholtz-like instabilities discussed above. We study the effect of varying the canopy parameters on this instability as well. A model based on linear stability analysis is proposed to capture some of these effects.
The paper is organised as follows. The numerical methods used and the canopy parameters are discussed in §2. The results from the DNSs, 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.

Methodology
We conduct direct numerical simulations of symmetric channels with rigid canopy elements on both walls. The streamwise, wall-normal and spanwise coordinates are x, y and z, and the associated velocities u, v and w. The wall-normal origin, y = 0, is defined at the tip plane of the canopies protruding from the bottom wall. The channel height, 2δ, 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δ 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πδ in the streamwise direction and πδ in the spanwise direction, with δ = 1. The domain-to-canopy height ratio for most cases considered here is (δ + h)/h ≈ 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. (2017) and MacDonald et al. (2018), 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 δ/s 10, and for the canopy with the largest element spacing is δ/s ≈ 4. The flow is incompressible with the density set to one. The simulations are run at a constant flow rate, with the viscosity adjusted to obtain a friction Reynolds number Re τ = u τ δ/ν ≈ 185, where u τ is the friction velocity calculated at the canopy tips. The simulation parameters are given in table 1 for reference.
The numerical method used to solve the three-dimensional Navier-Stokes equations is adapted from Fairhall & García-Mayoral (2018). 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 ∆y + min ≈ 0.33 at the canopy-tip plane, stretching to ∆y + max ≈ 3.3 at the channel centre. The grid within the canopies preserves the resolution of ∆y + min ≈ 0.33 near the canopy-tip plane, and for the tallest canopies considered stretches to ∆y + max ≈ 4 at the base of the canopy, where the flow is quiescent. 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 2011;Fairhall & García-Mayoral 2018;Abderrahaman-Elena et al. 2019). In the central block, the resolutions in the streamwise and spanwise directions are ∆x + ≈ 6 and ∆z + ≈ 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 fineblock 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 1991) where I is the identity matrix and L, D and G are the Laplacian, divergence and gradient operators respectively. N is the advective term which is dealiased using the 2/3-rule (Canuto et al. 2012). The Runge-Kutta coefficients, α k , β k , γ k and ζ k , for each substep, k, are taken from from Le & Moin (1991). The time step is ∆t. The canopy elements are represented using an immersed-boundary method adapted from García-Mayoral & Jiménez (2011). 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 topview 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 + ≈ 96, and element spacings ranging from s + ≈ 10 to 48. The second group, marked by the prefix 'H', consists of canopies with a fixed element spacing, s + ≈ 16, and element heights ranging from h + ≈ 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 ≈ 4, and w/s = 2/9. The heights for the canopies of G range from h + ≈ 10 to 100, with the element spacings varying in proportion to their height. Two additional simulations, H32 180 and H32 400 , are conducted to check the dependence of the results on the friction Reynolds number. The canopy geometries for both these simulations have s + ≈ 16, h + ≈ 32 and w/s = 1/2, with friction Reynolds numbers Re τ ≈ 180 and 400. We also conduct 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 rms 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 underpredict 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 H32 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 H32 180 and H32 400 to avoid grid related discrepancies in the comparison of their results.

Reynolds number effect
To analyse the influence of the Reynolds number in our subsequent DNSs, we compare the results of cases H32 180 and H32 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 + ≈ 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. 1999). Scaling in friction units over conventional rough surfaces has also been noted by Chan et al. (2015). Beyond y + 10, we observe that the magnitude of the peaks in the fluctuations and the Reynolds shear stresses are larger for case H32 400 compared to case H32 180 . The increase in magnitude of the near-wall peaks in the velocity fluctuations at friction Reynolds numbers larger than Re τ ≈ 180 is consistent with that observed in smooth-wall flows (Moser et al. 1999;Sillero et al. 2013), also included in figure 3 for reference. Further away from the canopy tips, at y + > 50, the rms 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 rms fluctuations being similar for these simulations, the distribution of energy in different scales is also similar. This is illustrated by the pre-multiplied spectral energy densities at y + ≈ 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 1999). The results of H32 180 and H32 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 ≈ 180 should also be relevant for higher Reynolds number flows.

Effect of canopy parameters on the surrounding turbulence
In this section, we discuss the effect of changing the height and spacing of the canopy elements on the velocity fluctuations within and above the canopies. 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 et al. 1991;Flack et al. 2007; Abderrahaman-Elena et al. 2019). Jiménez (2004) 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. (2017) for high-aspect-ratio prismatic roughness and by MacDonald et al. (2018) for spanwise aligned grooves with large spacings, and has been used to model cuboidal roughness by Yang et al. (2016). 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 (1972) where u is the full velocity, U is the mean velocity obtained by averaging the flow in time and space, and u is the full space-and time-fluctuating signal. The latter is decomposed into the element-induced velocity, u, which is obtained by averaging the flow in time alone, and the incoherent, background-turbulence fluctuating velocity u . We observe that the element-induced velocity 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 element spacing is the relevant lengthscale for the overlying flow. The element-induced fluctuations decay below the canopy-tip plane as well. The elementinduced fluctuations within become zero well above the canopy base for the canopies G10 and S10, that is, the canopies with the smallest element spacings. For the other canopies, the fluctuations only become zero at the wall at the canopy base. The intensity . Root-mean-square velocity fluctuations of the element-induced flow. The lines 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. of the velocity fluctuations within the canopy also depends on the element spacing. 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 figures 6(b, e, h). For the canopies of families S and G, the intensity of the velocity fluctuations within the canopy increases with element spacing, when scaled using either the friction velocity or the channel bulk velocity. This suggests that the intensity of the element-induced flow within the canopy also depends essentially on the canopy spacing.
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. Abderrahaman-Elena et al. (2019) observed a similar effect over conventional cubical roughness, where the element-induced fluctuations only extended to y ≈ h, but the effect of the roughness on the overlying flow extended to y ≈ 3h above them. At heights of y/s > 2-3 above the canopy-tip plane, the full rms 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 pre-multiplied spectral energy densities of the canopy and smoothwall 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 + 90. This corresponds to a height of about 2s for case S48, the canopy with the largest spacing. Although not shown here, the pre-multiplied spectral energy densities of the canopies of families H and G collapse with the smooth-wall spectra for y/s 3 as well.
A. Sharma and R. García-Mayoral Turbulent flows over dense filament canopies 1  Figure 7. Rms velocity fluctuations within and above the canopies. The lines 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. The black lines represent the smooth-wall case, SC.
kxkzEuu kxkzEvv kxkzEww kxkzEuv Figure 8. Pre-multiplied spectral energy densities at y + ≈ 90, with line contours from red to blue representing cases S10 to S48, normalised by their respective uτ . 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.
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 smoothwall like (Nikuradse 1933;Raupach et al. 1991). Typically, roughness elements with a characteristic size of a few wall-units, h + 5, fall into the hydraulically-smooth category The lines from red to blue 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. (Raupach et al. 1991;Jiménez 2004;Flack et al. 2007). Of the canopies studied here, we observe that the flow over the canopy of case G10, which has an element spacing of s + ≈ 2.6, is essentially smooth-wall like. This is evidenced by the collapse of the rms 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 smoothwall-like behaviour. 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 velocity fluctuations, 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 , typical of smooth-wall flows, is severely reduced. We also observe an increase in the Reynolds shear stresses above the canopy tip plane with increasing element spacing, shown in figures 9(d, f ), with an associated increase in the drag. The drag increase due to 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 1954). This shift can be observed for the canopies of families S and G in figures 9(a, c). The decrease in u and increase in v , w and drag with increasing element size is commonly reported over conventional rough surfaces (Ligrani & Moffat 1986;Orlandi & Leonardi 2006;Abderrahaman-Elena et al. 2019). Several authors have attributed the changes in the velocity fluctuations to the roughness elements modifying the near-wall turbulence cycle (Jiménez 2004;Flores & Jiménez 2006;Flack et al. 2007;Abderrahaman-Elena et al. 2019). The associated structures, such as streaks and quasi-streamwise vortices, are predominantly streamwise coherent (Kline et al. 1967;Jiménez & Pinelli 1999). A breakdown of the streamwise coherent structures in the flow with increasing canopy element spacing can also be observed in the present simulations. This is illustrated in the instantaneous realisations of the wall-normal velocity for the canopies of family G, portrayed in figure 5. Below the canopy-tip plane, increasing the element spacing results in an increase in all the components of the velocity fluctuations. The wall-parallel fluctuations, u and w , drop rapidly below the canopy-tip plane, and their magnitude reaches a plateau in the middle of the canopy, before they drop again near the canopy base to meet the no-slip condition. The height over which the fluctuations decay within the canopy and the magnitude of the fluctuations in the middle of the canopy appear to correlate with the element spacing. Note that this plateau in u and w 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 ≈ 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 fluctuations account for less than 10% of the magnitude of the total velocity fluctuations within the canopy.
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 rms velocity fluctuations for these cases. Above the canopy tips, we observe a decrease in u and an increase in v and w with an increase in the canopy height, similar to the effect of increasing the element spacing, as shown in figures 7(b, e, h). Within the canopy, u and w for all cases collapse to a single curve, only departing from it to meet the no-slip condition at the canopy base. The magnitude of the wall-normal fluctuations within the canopy, however, increases with the canopy height up to h/s ≈ 6. For h/s 6, the effect of the canopy height on the flow, within and above the canopy, saturates. This is 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 of cases H96 and H128 with h/s ≈ 6 and 8, respectively, are essentially the same.
The variations observed in the velocity fluctuations for the fixed-element-spacing simulations may, in part, result from the growth of a Kelvin-Helmholtz-like instability, typically reported in dense canopy flows (Finnigan 2000;Nepf 2012). In order to assess the presence of this instability in the flow, we compare the pre-multiplied spectral energy densities of the wall-normal velocity at y + ≈ 15 in figures 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, λ + z > 100, for a narrow range of streamwise wavelengths, λ + x ∈ (150-250). This range of streamwise wavelengths remains roughly constant for increasing canopy heights. Such a footprint in the spectral energy densities has been previously associated with the presence of spanwise-coherent, Kelvin-Helmholtz-like instabilities over riblets . Above a short canopy, like that of H16, the impermeability condition at the base of the canopy would inhibit the instability by blocking the wall-normal fluctuations (Huerre 1983;Healey 2009). Increasing the canopy height weakens this effect, leading to a stronger footprint of the instability, as observed in figures 10(f -j). 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  figure 9(e). The shift in the mean-velocity profile for these canopies, which reflects the increase in drag, can be observed in figure 9(b). The above discussion suggests that the growth of the Kelvin-Helmholtz-like instability is, in large part, responsible for the changes observed in the velocity fluctuations, Reynolds shear stresses and friction drag for the fixed-spacing canopies of family H, and that the effect of the element-induced flow is secondary. The growth of the Kelvin-Helmholtz-like instability near the canopy-tip plane 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 wallnormal spectral energy densities of the flow within the canopies, portrayed at y + ≈ −10 for all the canopies of family H, in figures 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 + ≈ −10 for this case also shows some overlapping regions 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 correspond to y + ≈ 1, which is as low as possible while yielding a non-negligible energy, since no direct comparison Figure 11. Pre-multiplied spectral energy densities of the wall-normal velocity, kxkzEvv, for (a-e) cases H16 to H128 at a height of y + ≈ −10; and (f -j) cases S10 to S48 at a height of y + ≈ −40. The contours are normalised by the rms values of their respective cases. The shaded contours are of the smooth-wall case, SC at a height of y + ≈ 1, for reference. The contours are in increments of 0.075 for all the cases.
with y + ≈ −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. 2019). 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 larger. 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 + ≈ −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 + ≈ −40 are portrayed in figures 11(f -j) and 12, respectively. It is also worth noting that even in canopies with small element spacings, such as that of case S10, the wall-parallel velocity fluctuations decay rapidly below the canopy-tip plane, but the wall-normal velocity fluctuations decay much slower, as shown in figure 7. This is also illustrated by the velocity fluctuations of the canopies of family H, where the wall-normal velocity fluctuations within the canopy require larger canopy heights to asymptote compared to the wall-parallel fluctuations. The presence of wall-normal fluctuations deep within the canopy are a reflection of the canopy layout obstructing the wall-normal flow less 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. Within a canopy with a large height, the only mechanism to inhibit the velocity fluctuations away from the canopy base is through the effect of the canopy drag. As the canopy geometries studied here exert more drag on the wall-parallel flow than the wall-normal flow, u and w decay faster than v within the canopy. The formation of the Kelvin-Helmholtz-like instability over canopies is also affected by the spacing between elements. To study this effect, we compare the pre-multiplied spectral energy densities of the wall-normal velocities for the canopies of family S, which have a constant height and different element spacings. For the canopy with the smallest spacing, S10, the energetic wavelengths in the flow at y + ≈ 15 are similar to those in smooth-wall flows, and the footprint of the instability in the flow is weak, as shown in figure 10(a). As the element spacing is increased, there is an increase in the energy in wavelengths associated with the Kelvin-Helmholtz-like instability. This suggests that this instability is inhibited in very dense canopies, possibly due to the large drag they exert on the velocity fluctuations Sharma et al. 2017). 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 large streamwise wavelengths λ + x 700. We also observe that, while in the canopies of family H the streamwise wavelength of the instability was roughly constant independently of the canopy height, the increase in the element spacing results in an increase in the streamwise wavelength of the instability from λ + x ∼ 140 for case S10 to λ + x ∼ 280 for case S48. 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 (Poggi et al. 2004;Pietri et al. 2009;Huang et al. 2009;Sharma & García-Mayoral 2019).
Finally, let us focus on the self-similar canopy geometries of family G. The effect of increasing the size of the canopy elements produces similar effects on the pre-multiplied spectral energy densities as for the canopies of family S, discussed in the previous paragraph. For the densest canopy, G10, the spectral energy densities at y + ≈ 15 collapse with those over a smooth wall, as shown in figures 13(a-d). As the size of the canopy is increased, we observe a stronger footprint of the Kelvin-Helmholz-like instability in the energy densities, and there is also a progressive increase in the streamwise wavelengths associated with the instability. We also 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 apparent in the energy densities of the other velocity fluctuations and the Reynolds shear stresses. 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 large spanwise wavelengths compared to smoothwall flows. We can also observe these effects qualitatively in the instantaneous realisations of the wall-normal velocity, shown in figure 5. For case G10, the structures observed in the instantaneous flow field are predominantly streamwise-coherent, as for smooth walls (Kline et al. 1967;Jiménez & Pinelli 1999). As the canopy size increases, we observe a gradual breakup of this streamwise coherence and an increase in the presence of spanwisecoherent structures.
The results discussed in this section suggest that the growth of the Kelvin-Helmholtzlike 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. The streamwise wavelength of the Kelvin-Helmholtz-like instability is determined by the shear length, L s = U/(dU/dy), calculated at the canopy-tip plane (Raupach et al. 1996;Nepf 2012). Nepf et al. (2007) showed that the shear length L s in canopy flows, in turn, is determined by the mean streamwise canopy drag coefficient. Intuitively, in tall, dense canopies, we would expect this drag coefficient to be a function of the element spacing. Therefore, the canopies of family H, which have a constant element spacing, will 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 would decrease the canopy drag coefficient, thereby resulting in the larger instability wavelengths observed. The effect of the canopy spacing on the drag coefficients and the instability wavelengths will be discussed further in §4. With regard to the effect of the canopy parameters on the intensity of the instability,  postulated that the growth of the instability over dense canopies was governed by two competing effects, the shear at the canopy tips and the canopy drag. Canopies with closely packed elements would result in larger shear and a stronger inflection point at the canopy-tip plane, which would enhance the instability, but at the same time, a denser canopy would also produce a larger drag, which would inhibit the instability (Sharma et al. 2017). This effect is observed for the canopies of families S and G, where increasing the element spacing, or reducing the canopy drag, results in a stronger footprint of the instability in the flow.

Effect of Re τ 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 friction Reynolds numbers. In this section, we expand discussion in §2.1 to discuss the effect of changing the friction Reynolds number on the Kelvin-Helmholtz-like instability over dense canopy flows. We have noted previously that the presence of this instability in the flow leaves a clear footprint in the wall-normal spectral energy densities, in the form of additional energy in large spanwise wavelengths for a distinct range of streamwise wavelengths. It can be observed in figure 4(b) that the footprint of this instability in the simulations of H32 180 and H32 400 are essentially the same, and the streamwise wavelength of the instability for both cases is λ + x ≈ 150. This suggests that the instability, like the turbulent fluctuations, also scales in friction units. The scaling of the Kelvin-Helmholtz-like instability in friction units has also been noted over riblets (García-Mayoral & Jiménez 2012) and some permeable substrates (Gómez-de-Segura & García-Mayoral 2019). At the same time, the wavelength and amplification of the instability are governed by the shear at the canopy tips (Raupach et al. 1996;Ghisalberti & Nepf 2002). For very dense canopies such as those considered here, the shear at the canopy-tip plane, y = 0, also determines the friction velocity, which explains the common scale for the background turbulence and the instability.

Linear analysis of Kelvin-Helmholtz-like instabilities
The results from DNS discussed in §3 show that the flow in the region near the canopytip 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 twodimensional, mean-flow linear stability analysis, even in turbulent flows (Jimenez et al. 2001;García-Mayoral & Jiménez 2011;Gómez-de-Segura & García-Mayoral 2019). In this section, we discuss the methodology and results from a temporal, mean-flow stability analysis on the velocity profiles obtained from the DNSs. As the Kelvin-Helmholtz-like instability is an inviscid phenomenon, several of the aforementioned studies 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. 2001;Gómez-de-Segura & García-Mayoral 2019). 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 a stability analysis performed on modelled a priori velocity profiles which would not require any information from the DNSs.
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 where C i is the effective canopy drag coefficient in each i th direction, assumed to be homogeneous over the entire canopy, as in Singh et al. (2016). 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 (Tanino & Nepf 2008).
In the centre of the canopy elements, 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 The colours from red to blue represent cases S10 to S48, H16 to H128 and G10 to G100.
1856), and has been used by  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 dP /dx obtained from the DNSs into equation (4.3). From dimensional arguments, equation (4.3) predicts that the drag coefficient would scale as C x ∼ ν/s 2 . This scaling is demonstrated in figure 14(a), which suggests that equation (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 (∂ 2 x + ∂ 2 z )v = dP/dy. The wall-normal drag coefficient is then obtained as 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 ≈ 2-3, which shows that the streamwise flow is more obstructed than the wall-normal flow for the layouts considered.
Linearising the equations (4.1) and (4.2) around the mean flow, U (y), yields These equations are used to obtain a modified Orr-Sommerfeld equation (Drazin & Reid 2004;Singh et al. 2016;, (d-f ) mean profiles obtained from DNSs, with no drag on the perturbations; and (g-i) mean velocity profiles obtained using equation (4.10), with no drag on the perturbations. The lines 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. eigenvalue problem where the prime superscript denotes differentiation with respect to y, and D represents the operator d/dy. Equation (4.9) is solved to obtain ω for real values of α, subject to no-slip and impermeability boundary conditions at the top and bottom walls. The instability is amplified for positive values of the imaginary part of ω.
The most amplified wavelengths predicted by the stability analysis only match those observed in the DNSs for canopies with high values of δ/h. The growth rates for different perturbation wavelengths are portrayed in figures 15(a-c), and the wavelengths with the highest growth rates are summarised in table 2. The instability wavelengths predicted for cases H16, H32, G10, G20 and G40 show reasonable agreement with those observed in the DNSs. For canopies with larger heights, however, the analysis predicts wavelengths larger than those observed in the DNSs. For the fixed-spacing canopies of family H, the predicted instability wavelength also increases with increasing canopy height, whereas the DNSs 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  Table 2. Most amplified instability wavelengths observed in the DNSs 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. SAC0, most amplified wavelengths from stability analysis on DNS mean profiles without drag on fluctuations; SAMC0, on synthesised velocity profiles without drag on fluctuations; and SA, on DNS mean profiles with drag on fluctuations.
amplified wavelength, λ + x ≈ 385, portrayed in figure 16(a), show that it has a large wall-normal span, extending up to y + ≈ 120. Such an instability was also reported by Singh et al. (2016), who performed stability analyses similar to the one conducted here, except that the canopy was represented by a drag force varying with the square of the velocity. Singh et al. (2016) 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 δ/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 width of the channel. It is assumed in the present stability analysis that the drag coefficient experienced by the perturbations is the same as that experienced by the mean flow. It is possible, however, that scales much smaller than the mean, such as those of the instability, perceive a smaller drag coefficient, as reported in Sharma & García-Mayoral (2019) for sparse canopies, thus not exciting this secondary instability. Excluding the drag on the perturbations in the stability analysis yields better estimates for the instability wavelengths observed in the DNSs, as shown in figure  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 DNS 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. There may also be some distortion of the instability by the ambient turbulent fluctuations in the DNSs (Rogers & Moser 1994;Raupach et al. 1996). The stability analysis also indicates that the observed instability scales in friction units, as discussed in §2.1. Consistent with the observations in the DNSs, the linear stability analysis for cases H32 180 and H32 400 also predicts similar instability wavelengths and growth rates in friction units, as shown in figure 17.
The results obtained from the DNSs 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 shear-layer thickness within the canopy is commonly defined as L s = U/(dU/dy) at the canopy tips (Raupach et al. 1996;Finnigan 2000;Nepf 2012). However, this definition misses the contribution of the part of the shear layer above the canopy. García-Mayoral & Jiménez (2011) studied the formation of Kelvin-Helmholtz-like instabilities over riblets, noting that the shear-layer thickness above was given by the height at which the vorticity gradient, d 2 U/dy 2 , concentrated. In smooth-wall flows, this height is roughly y + c ≈ 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 , taking y + c = 5, as shown in figure 18(a). This suggests  Figure 18. (a) Instability wavelength, λ + x , obtained from the linear stability analysis versus the total shear length, L + s + y + c ; (b) shear length, L + s , versus the drag lengthscale; (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.
that the shear-layer semi-thickness above the canopies is roughly constant for most of the geometries considered here, and remains close to the smooth-wall value. The only notable deviation is for the sparsest canopy studied, S48. For canopies with large element spacings, we observe that the peak in d 2 U/dy 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. We observe that the height of the shear layer within the canopy, L s , is set by the mean canopy drag coefficient, L s ∝ ν/C x , which in turn depends mainly on the element spacing, 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 results and the stability analysis.

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, the turbulent stresses do not penetrate within , 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. 2001;García-Mayoral & Jiménez 2011;Gómez-de-Segura et al. 2018;Gómez-de-Segura & García-Mayoral 2019). The equation for the mean velocity can then be written where C x (y) is the average streamwise canopy drag coefficient, which is assumed constant within the canopy and zero outside, and ν T (y) is the height-dependent eddy viscosity proposed by Cess (1958) 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 equation (4.3) and portrayed in figure 14(a). These have been obtained using the data from the direct numerical simulations, but can be obtained from Stokesflow simulations which are less computationally intensive, similar to how we obtained the wall-normal canopy drag coefficient, C y . The most amplified wavelengths predicted by the stability analysis conducted on these modelled velocity profiles, with no drag on the fluctuations, are in reasonable agreement with those conducted on profiles obtained from the DNSs. The growth rates predicted are portrayed in figures 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. 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 DNSs, 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, which fails as the element spacing is increased. As discussed previously, the wavelength of the instability is set by the shear length. The shear length within the canopy, L s , is set by the canopy drag coefficient, C x . As this drag coefficient is the same both from the DNSs and for the modelled velocity profiles, we expect L s to be similar as well. The shear length above the canopy, however, may differ, as the profiles from DNS 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 the dense canopies considered in this work, turbulence is essentially precluded from penetrating into the canopy and the shear length above does not seem to vary significantly from its smooth-wall value.

Conclusions
In the present work, we have examined the effect of canopy geometry 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 ranging from s + ≈ 2.6 to 48. We have observed that the height of the roughness sublayer in such canopies is set by the element spacing, and not their height. The intensity of the velocity fluctuations within the canopy is also determined by the element spacing. The element-induced fluctuations are essentially the same for canopies with the same element spacing regardless of their height, and increase in magnitude with the spacing. As the element spacing increases, so do the peak intensity of the wall-normal and spanwise velocity fluctuations and the Reynolds shear stresses above the canopy. The peak in the streamwise velocity fluctuations, in turn, decreases with increasing element spacing. These changes can be attributed to the modification of the near-wall turbulence through its interaction with the canopy elements, and to the growth of a Kelvin-Helmholtz-like instability at the top of the canopy. At a fixed element spacing, the effect of increasing the height of the canopy is found to asymptote at a height-to-spacing ratio h/s ≈ 6. While increasing the canopy height at a fixed spacing does not have a significant effect on the element-induced flow, it promotes the formation of the Kelvin-Helmholtz-like instability, which in turn affects the turbulence fluctuations within and above the canopy. We also find that the velocity fluctuations deep within the canopy result largely from spanwisecoherent scales, suggesting that they originate from the overlying instability as well.
We have also presented results from linear stability analysis conducted using the mean velocity profiles obtained from the DNSs. This analysis is able to capture the approximate wavelength of the instability observed in the DNSs 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 show that the shear length within the canopy, which determines the instability wavelength, increases linearly with the element spacing, and therefore so does the instability wavelength. We have proposed a simplified model to capture the most amplified 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.
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. , 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.