The role of the in-plane solidity on canopy flows

Abstract Turbulent open channel flows developing above submerged canopies made of slender cylinders mounted perpendicular to the channel bed are known to be largely governed by the solidity parameter $\lambda =dh/\Delta S^2$ ($d$ and $h$ being the diameter and height of the filament, and $\Delta S$ the average spacing between filaments). When the filaments are sufficiently slender, the ratio between the height of the stems and the spacing sets the hydrodynamic regime developing inside and outside the canopy. This ratio also establishes the conditions leading to the transition from a dense to a sparse canopy flow regime (Nepf, Annu. Rev. Fluid Mech., vol. 44, 2012, pp. 123–142). In a previous, companion numerical investigation, Monti et al. (J. Fluid Mech., vol. 891, 2020, A9) used large eddy simulation (LES) to study the influence of the canopy height on the onset of the different regimes without modifying the average spacing $\Delta S$ between the stems. In that LES study, we were looking at the complementary situation in which the height of the stems is constant while the filaments’ number density of the canopy is changed. It was found that for low values of $\lambda$ (i.e. sparse or moderately dense canopies: $\lambda \lessapprox 0.26$), the flows sharing the value of the solidity obtained by either varying $h$ or $\Delta S$ are very similar. Differently, for higher values of $\lambda$ (i.e. in denser canopies), the effects of $h$ and $\Delta S$ start to diverge although sharing the same nominal value of $\lambda$. In this paper, we analyse the different physical mechanisms that come into play for dense configurations obtained by varying either $\Delta S$ or $h$. In particular, we focus on the most relevant length scales and carry out a detailed analysis of the flows using a triple decomposition approach. We show that the inner region of dense canopy flows, characterised by tall stems, is dominated by wall-normal sweeps delivering high momentum in the wall vicinity. Here, the impenetrability condition of the bed redistributes the available momentum in the wall-parallel directions re-energising an otherwise stagnating flow. Differently, in densely packed canopies, the penetration of the outer jet and the momentum transfer from the external flow are limited by the decreasing value of the wall-parallel permeabilities leading to different behaviours, including a reduction of the total drag offered by the canopy.


Introduction
Canopy flows play a crucial role in various natural and technological processes, making them a subject of significant interest.They have been studied in the context of geophysical flows over crops or forest vegetation and aquatic canopies.One important parameter characterising canopy flows is the ratio between the depth of the canopy (average vegetation height, h) and the thickness of the boundary layer (δ) that forms above it.For fully submerged canopies in aquatic environments, the canopy length scale is represented by the open channel or river depth H.In this study, we focus on canopies with a fixed ratio of H/h = 10, which falls under the category of 'fully submerged canopy' according to the classification proposed by Nepf (2012).
In the context of this class of canopies, Nepf (2012) introduced a further geometric classification based on the concept of solidity (λ).The latter is defined as the ratio of the frontal projected area of the canopy to its base area and, thus, provides general information about the arrangement and spacing of the canopy elements.In this study, we consider an idealised canopy composed of cylindrical filaments, as shown in figure 1.In this case, the parameter λ can be expressed as where a(y) = d(y)/ S 2 is the frontal area per canopy volume, d is the diameter of the cylinders (the filaments), h is the canopy height and S is the mean spacing between filaments.
Previous research by Nepf (2012) and Poggi et al. (2004) has suggested that the value of λ can be used as an indicator to predict the flow regime developing in canopy flows.For λ < 0.1, the flow develops into a sparse regime where the drag offered by the canopy is small compared to the one induced by the bed.In contrast, for λ > 0.1, the drag offered by the canopy becomes dominant.The abrupt discontinuity in the drag profile at the canopy tip is responsible for the appearance of an inflection point in the mean velocity distribution and may induce a Kelvin-Helmholtz (KH) instability, generating large-scale spanwise coherent motions (Raupach, Finnigan & Brunet 1996;Finnigan 2000;Ghisalberti & Nepf 2002).
Several conceptual models have been proposed to explain the emergence of these spanwise coherent structures above dense canopies, relating them to the roll-up of Kelvin-Helmholtz waves induced by the inflected mean velocity profile.Raupach et al. (1996) observed that the streamwise wavelength of the most amplified Kelvin-Helmholtz mode Λ x is preserved downstream of the canopy in fully developed conditions.Another length scale, known as the shear length scale (L s ), was introduced to characterise the vorticity layer thickness above the canopy and to introduce an analogy with mixing layer flows.The shear length scale, defined as the ratio of the canopy tip velocity to the wall-normal gradient of the mean velocity U(y) (Raupach et al. 1996), Figure 1.Geometrical parameters that characterise an idealised canopy composed of cylindrical filaments (Nepf 2012).The velocity profile (given in blue) corresponds to the one of a typical dense canopy flow regime (Monti, Omidyeganeh & Pinelli 2019).
is proportional to the thickness of the vorticity layer δ ω .Raupach et al. (1996) hypothesised that in dense regimes, the Λ x /L s ratio falls between 7 and 10.Multiple experimental datasets substantiated this range, further narrowing down Λ x to an estimate of 8.1L s .Ghisalberti & Nepf (2004) reported consistent thickness of the vorticity layer in the streamwise direction above the canopy, validating the fully developed mixing layer analogy.However, the 'roller'-shaped structures above the canopy exhibited an undefined vorticity pattern, complicating their direct detection (Finnigan, Shaw & Patton 2009).
Several studies using linear stability analysis have shed light on the temporal evolution of instabilities leading to canopy-top structures (White & Nepf 2007;Luminari, Airiau & Bottaro 2016;Zampogna et al. 2016).Another study reported an intensification of the KH instability signature with increasing mean filament spacing ( S), assuming a linear relationship with the shear length scale (L s ) (Sharma & García-Mayoral 2020b).Monti et al. (2019) have recently observed how the outer, quasi-streamwise vortices cause the loss of the spanwise coherent nature of the vorticity rollers.
In recent work, Monti et al. (2020) identified novel coherent motions within dense canopies, potentially originating from an instability at the inner inflection point of the mean velocity profile.These coherent motions, combined with the wall-normal penetration of outer structures, significantly impact the modulation of velocity components within the canopy.
Canopy flows are typically analysed using a zonal approach, as suggested by various authors (Belcher, Jerram & Hunt 2003;Poggi et al. 2004;Okamoto & Nezu 2009;Monti et al. 2020).This approach categorises each wall-parallel slice of the flow according to its specific behaviour, offering a comprehensive understanding of flow dynamics within the canopy.
The question of the conditions that trigger different flow behaviours across canopy layers and their modulation by respective canopy regimes, particularly in the transition regime, remains open (Nicholas, Omidyeganeh & Pinelli 2022).To tackle this and other questions involving complex dynamics, our study uses stem-resolved simulations.These high-fidelity simulations, capable of capturing the intricacies of flow structures within the canopy, present a valuable alternative to traditional experimental methods and numerical simulations based on drag models.
Although earlier drag model-based simulations provided reasonable results in capturing the flow behaviour near the canopy, they fell short in providing a detailed characterisation of the flow within the canopy layer and at the interface with the outer flow.Recent advances in immersed boundary formulations have facilitated the use of stem-resolved simulations, which have proven successful in capturing the dynamics within the canopy layer (Monti et al. 2019;Sharma & García-Mayoral 2020a).Recently, using this methodology, ) computed at the virtual origin (to be introduced later).Additionally, y + w and y + c are the grid spacings, in viscous units, within the region next to the impermeable wall.Monti et al. (2020) were able to explore the inception of different flow regimes when changing the canopy's frontal solidity.In this work, we will be focusing on the complementary problem that concerns the impact of the in-plane solidity on the genesis of different flow regimes.The manuscript is organised into several sections: § 2 describes the geometrical configurations and numerical techniques used for the simulations, while § 3 analyses five different canopy configurations (see table 1) using both statistical characterization and flow field instantaneous realisations.A few conclusions are presented in § 4.

Geometry set-up and numerical procedure
We performed high-fidelity simulations of turbulent flows over rigid canopies using our in-house developed incompressible Navier-Stokes solver called SUSA (Omidyeganeh & Piomelli 2013a;Monti et al. 2019).The simulations were conducted using a resolved LES approach, where the governing equations are obtained by filtering the velocity and pressure fluctuations below a spatial filter with a typical size in the inertial sub-range of turbulence (Moser, Haering & Yalla 2021).
In the Cartesian framework, the streamwise, wall-normal and spanwise directions are represented by x, y and z (or sometimes denoted as x 1 , x 2 and x 3 ), respectively.The velocity components along these directions are denoted by u, v and w (or u 1 , u 2 and u 3 ), and the pressure is denoted by p.
The governing LES equations in dimensionless form are given by Here, Re b = U b H/ν represents the bulk Reynolds number based on the bulk velocity (U b ), open channel height (H) and the kinematic viscosity (ν).The sub-grid Reynolds stress tensor, τ i,j = u i u j − u i u j (Leonard 1975), is modelled using the integral length scale approximation (ILSA) method (Piomelli, Rouhi & Geurts 2015).In this approach, the length scale and model constant required for the closure problem are computed locally, eliminating dependence on the grid (Rouhi, Piomelli & Geurts 2016).
The spatial discretisation of the incompressible LES equations is based on a second-order accurate, cell-centred finite-volume formulation.The pressure p and velocity components u i are co-located at the centre of the cell to avoid spurious pressure modes, employing a deferred correction approach (Rhie & Chow 1983).
The time advancement of the governing equations is performed using a second-order, semi-implicit fractional step procedure (Kim & Moin 1985;Chorin 1968).The nonlinear and wall-parallel diffusive terms are treated using a second-order accurate, explicit Adam-Bashforth scheme, while the wall-normal diffusion is handled with a second-order accurate, implicit Crank-Nicolson scheme.
Regarding the canopy treatment, in contrast to other studies (Finnigan et al. 2009;Bailey & Stoll 2016), we resolve the filamentous layer at the stem level by imposing a zero-velocity condition along each slender, rigid cylindrical rod composing the canopy.We employ the immersed boundary method (IBM) (Pinelli et al. 2010) to enforce the zero-velocity condition on each individual filament.
For additional information on the code, its parallelisation and the validation conducted for other flow configurations, refer to the following publications: Omidyeganeh & Piomelli (2011, 2013a,b); Rosti, Omidyeganeh & Pinelli (2016).The role of sub-grid scale stresses, and the behaviour of the LES filter operator in the canopy and outer layer are discussed by Monti et al. (2020).For detailed descriptions of the IBM, including its properties and implementation for various flow configurations, consult Pinelli et al. (2010) and Favier, Revell & Pinelli (2014).The IBM assessment in the context of canopy flows, including calibration and extensive validation, is presented by Monti et al. (2019).In the same reference, the reader can also find a detailed description of the filament resolution: note that the method is able to capture the hydrodynamic diameter of the filaments but does not resolve their surface.
The five cases considered in this study share a common computational domain as the one considered in the study by Monti et al. (2020), we consider a half-channel with streamwise, wall-normal and spanwise directions lengths set to L x /H = 2π, L y /H = 1 and L z /H = 3π/2, respectively (see figure 2).The appropriateness of the chosen domain size was verified by Monti et al. (2019), who demonstrated that the large-scale turbulent structures are adequately captured within the computational box (i.e.velocity correlations decrease to small values at half the domain size in the streamwise and spanwise directions).
For all simulations, the origin in the wall-normal direction is set at y/H = 0, representing the location of the impermeable bottom wall, referred to as the canopy bed.With this set-up, the filaments form a canopy with a height of h/H.
To model an open channel flow, we enforce a zero-velocity condition at y/H = 0 and a free-slip condition (∂ y u = ∂ y w = 0, v = 0) at the top surface (y/H = 1).In the homogeneous x and z directions, periodic boundary conditions are applied.This choice is justified by the experimental observations of Ghisalberti & Nepf (2004) who observed a constant thickness layer developing above the canopy tip.
At each time step, a uniform pressure gradient is computed and imposed in the streamwise direction to enforce a constant, time-independent volumetric flow rate.The latter results in a bulk Reynolds number of Re b = U b H/ν = 6000.This choice allows for direct comparison with the results of Monti et al. (2020) and closely matches the experimental conditions reported in the works of Shimizu et al. (1991) and Ghisalberti & Nepf (2004).
Following similar previous numerical simulations (Thakkar, Busse & Sandham 2018;Monti et al. 2020), the distribution of the filaments is carried out using a tiled approach: the bottom impermeable wall is divided into non-overlapping square tiles with a uniform edge length of S (see figure 1).Within each tile, a single filament is randomly positioned, using a uniform distribution.This approach prevents the formation of filament clusters or repetitive patterns that could lead to channelling effects within the canopy layer (Stoesser et al. 2009;Sharma & García-Mayoral 2020a).
In the present study (summarised in table 1), we change the solidity λ = (hd)/ S 2 by altering only the tile size ( S), while keeping the canopy height and diameter constant (i.e.h/H = 0.1 and d/H = 0.024).Therefore, our investigation complements the work by Monti et al. (2020), where the canopy height was adjusted for a fixed average filament spacing.The chosen tile sizes in our study are selected to encompass a range of canopy flows, spanning from a sparse to a dense regime, as previously described by Nepf (2012).
Regarding the spatial discretisation, we employ a Cartesian mesh with a uniform distribution of nodes in the wall-parallel directions and a stretched arrangement along the wall-normal coordinate.This arrangement clusters grid points near the canopy bed and in the vicinity of the tip region.To ensure consistency across all configurations, a single computational mesh is determined based on the estimated boundary layer thickness (δ ≈ H − h).The mesh consists of N x = 576 grid points along the x direction, N y = 230 along the y direction and N z = 432 along the z direction, as summarised in table 1.The grid resolutions in wall units, based on the friction velocity computed using the total stress at the virtual origin (the origin observed by the outer flow), are presented in table 1 (where subscript 'o' indicates non-dimensional parameters based on the frictional velocity at the virtual origin).
The resolution in the wall-parallel directions, expressed in outer flow units (i.e. x + o = x × u τ,o /ν), satisfies the condition x + o = z + o ≤ 9 for all configurations.In the wall-normal direction, a stretching function is employed to maintain a resolution of y + o < 0.4 in the localised layers near the impermeable wall and encompassing the canopy tip.
The chosen grid resolution adheres to standard values used in DNS studies of typical wall-bounded turbulent flows (Kim, Moin & Moser 1987;Lee & Moser 2015).If the canopy bed friction velocity is used to determine the wall units, the mesh provides a resolution of x + i = z + i ≤ 3 and y + i ≤ 0.16, which also aligns with accepted DNS standards.
For further details on the numerical grid distribution, mesh convergence analysis and the suitability of the numerical scheme, refer to Monti et al. (2019).All the results presented in the subsequent sections are obtained from statistically converged simulations, with statistical data collected over a duration of 900tU b /H time units, approximately corresponding to 150 full washouts.Finally, regarding the comparison that will be presented in figure 24 ( § 3.4), which features an open channel bounded by a smooth wall at Re τ = 710, the results were acquired using the same code.No filaments were involved, and an identical domain was maintained with grid spacings x + = z + = 7.1 and y + (2) = 0.48 (where y(2) represents the coordinate of the first interior node in the wall-normal direction).

Results
In this section, we present the results obtained from simulations of five canopy flows defined by the parameters in table 1.The focus of the discussion is a direct comparison of flow statistics and structures between the cases studied here and those obtained by Monti et al. (2020).
Initially, we analyse the basic averaged quantities and their variations with respect to the in-plane solidity.We also assess the reliability of well-known non-dimensional parameters for identifying the canopy flow regime, particularly in relation to the shape of the mean velocity profile (Nepf 2012).
Subsequent subsections focus on higher-order statistics and the analysis of coherent structures that emerge or disappear with different canopy flow regimes.

Mean flow characteristics
We begin by examining the effect of in-plane solidity on the mean velocity profile.Distinct variations in the mean velocity distribution are evident for different solidity values λ, particularly within the canopy region (y/H ∈ [0, 0.1]).As the solidity increases, the mean velocity profiles exhibit increased retardation with a progressively enhanced concave shape.The sparse (SP) case demonstrates a minor slowdown, while the TR, MD, DE cases (i.e.cases with increasing solidity) exhibit mild retardation.In contrast, the densest (HD) configuration demonstrates a substantial decrease in mean velocity accompanied by significant curvature enhancement.
Comparing the present HD profile with the one obtained by Monti et al. (2020) at the same λ value (see figure 3b), clear differences in shape are evident.Specifically, the velocity gradient in the wall region is much milder in the present case, with a smoother change of curvature when moving away from the wall.The strong gradient in the wall region and the appearance of a kink in the velocity profile, just below the inner inflection point, reported by Monti et al. (2020), is absent in the present mean velocity profile.
The presence of the kink was attributed to the strong inrushes (jet-like) of wall-normal velocity (Banyassady & Piomelli 2015), deflected in the x and z directions due to continuity and wall impermeability.Although the two dense cases share the same solidity λ, the absence or presence of the kink suggests potential differences in the flow dynamics within the canopy region.
Furthermore, it is interesting to note that the mean velocity profile of the present HD case is similar to those reported in the context of turbulent wall flows bounded by porous surfaces (Breugem, Boersma & Uittenbogaard 2006;Rosti, Cortelezzi & Quadrio 2015;Kuwata & Suga 2017).
When considering a non sparse canopy flow regime (i.e.λ > 0.04), the mean velocity profile features three notable points within the filamentous layer (i.e. for y ≤ h as shown in figure 3a): two inflection points and the location of the virtual origin of the outer flow (Poggi et al. 2004;Nepf 2012;Monti et al. 2019).The latter can be interpreted as the location of the virtual wall seen by the outer turbulent boundary layer above the canopy.The concept of a shifted boundary layer is commonly proposed in the framework of boundary layers developing over textured surfaces (Jiménez 2004;Bottaro 2019).In our case, to determine the location of the virtual origin (y vo ), it is assumed that the outer mean velocity profile follows the canonical logarithmic law (Monti et al. 2020): where y out = y − y vo is the shifted wall-normal coordinate (i.e.y measured from the virtual origin), κ is the von Kármán constant and B is a shift accounting for the nature of the wall roughness (Jiménez 2004).In the present work, the friction velocity u τ,out in the above equation is computed from the total stress at the location of the virtual origin as the sum of the viscous stresses and the Reynolds shear stress: By ensuring the alignment of the von Kármán constant κ (here assumed to be κ = 0.41), the wall-normal position of the virtual origin for each canopy configuration can be determined if the mean velocity and total stress profile are known.The intercept value B in the logarithmic law is not an unknown, instead, it is a result obtained from solving (3.1) for y vo .Additional information regarding this process, as well as the sensitivity of the virtual origin's location to the chosen value of κ, can be found in the study by Monti et al. (2019).Figure 3(a) reveals the presence of two inflection points on each mean velocity profile: an outer inflection point located at the top of the canopy and an internal inflection point near the wall.The outer inflection point is generated by the abrupt discontinuity of the drag force exerted by the canopy at its tip.The internal inflection point results from the merging of the velocity profile close to the wall with the mean velocity distribution in the region closer to the canopy tip (Monti et al. 2020).
The position of the inner inflection point can be determined by finding the location where the second derivative of the mean velocity is zero.This can be accomplished by considering the mean momentum equation in the streamwise direction: 3), the first, second and third terms represent the mean streamwise pressure gradient, the viscous force and the Reynolds shear stress, respectively.The final term accounts for the overall mean drag imposed by the canopy filaments and, since the filaments are mounted normally to the bottom wall, its contribution is mainly due to the form drag of the cylinders, with a weak contribution from the viscous shear stress.
Figure 4 illustrates the mean locations of the notable points of the mean velocity profiles, namely the inner inflection point (red horizontal line), virtual origin (black horizontal line) and outer inflection point (blue horizontal line), as a function of the solidity λ for the five canopy configurations considered in this work.Note that when the virtual origin moves below the inner inflection point, no overlapping region can be defined any longer.
In the case of sparse canopies, characterised by S/h 1 (Poggi et al. 2004; Sharma & García-Mayoral 2020a), as λ is decreased, the virtual origin moves towards the canopy bed (i.e.y vo /H → 0) and the internal inflection point collapses onto the canopy tip.Conversely, as λ is increased, the location of the virtual origin starts moving away from the canopy bed and the inner inflection point moves towards it.The signed distance between these two points has a zero in the range 0.14 < λ < 0.25.The value of λ at which the two points collapse into a single location can be used as a criterion for predicting the onset of the dense canopy flow regime (Monti et al. 2020).
In dense configurations, the virtual origin is located in the upper half of the canopy (i.e.y vo /h > 0.5), while the inner inflection point is in the bottom half of the canopy (i.e.y vo /h < 0.5).An increase in the width of the overlap region suggests a decoupling between the inner and outer layers (Monti et al. 2020;Sharma & García-Mayoral 2020b).In the densest, HD case, the virtual origin and internal inflection points are located near the canopy tip and bed, respectively.In this situation, the outer flow behaves like a shear flow over a rough wall, while the innermost flow develops a local boundary layer near the bed. Figure 5(a,b) show the locations of the inner inflection point and virtual origin (the latter measured from the canopy tip) as a function of λ.These figures, which incorporate data from Monti et al. (2020), demonstrate that the parameter h − y vo and the location of the internal inflection point saturate to a constant value as λ increases (denser canopies).This behaviour can be better understood by considering a scenario where the filament spacing is kept constant (i.e.S/H = π/24) while systematically increasing the canopy height.In this case, the inner and outer layers gradually decouple until reaching conditions where the canopy height (or the solidity λ in this case) becomes irrelevant to the flow behaviour.
However, in the current framework where λ is increased by considering progressively denser canopies of the same height, no clear saturation is observed.The latter will eventually occur when extremely packed filaments are considered, i.e. almost solid wall as λ → ∞.
Further information on the locations of the notable points can be extracted by looking at their scaling with the canopy geometric parameters.Figure 6 shows the locations of the internal inflection points y in of the mean velocity profiles, scaled with S and h, plotted versus λ. Figure 6(a) demonstrates a collapse between the cases sharing the same λ value, regardless of whether it is obtained by changing S or h.As S approaches zero (i.e.increasing λ), the wall becomes solid and y in /h approaches zero as well, resulting in a linear decreasing trend.However, when h becomes very large, the behaviour of the inner region flow becomes independent of h (or λ) and y/h reaches an asymptotic value.
Figure 7 provides insights into the mean location of the virtual origin measured from the canopy tip as a function of λ.Unlike the thickness of the inner layer, which scales with S, the thickness of the outer layer scales with the canopy height h.The open symbols in figure 7(a) indicate that the outer thickness remains roughly the same as the filament average spacing, indicating that it is S that determines the size of the eddies able to penetrate the canopy (Monti et al. 2020).
However, when S is kept constant, the penetration is governed by h, and asymptotically saturates to unity as the canopy height increases towards the densest regime.Figure 7(b) shows that when (h − y vo ) is scaled with h, a collapse is observed irrespective of how the value of λ is obtained.This result is consistent with the interpretation of (h − y vo )/h as the fraction of the canopy that the outer eddies can penetrate.In summary, S determines the size of the eddies that can enter the canopy, while h determines the depth of penetration.

Notable points locations: heuristic model
We propose a simple geometric model to explain the variations in the position of the virtual origin of the mean velocity profile (Monti et al. 2020).Additionally, we discuss the influence of the mean filament spacing on the upwash/downwash events produced inside the canopy by outer turbulent eddies and their role in determining the transition between different canopy flow regimes (Nepf 2012).By keeping the canopy height h constant, the variations in S can be interpreted as modulations of the size of the high-pass filter acting on the outer flow structures.In particular, as illustrated in figure 8, it is the ratio ( S − d)/h that determines the size of the eddies capable of penetrating the canopy and their depth of penetration.For slender canopies (d/h 1), such as the cases under consideration in this work, the ratio S/h is as a sufficient indicator.
In the case of S/h > 1, when the width of the canopy voids exceeds the characteristic eddy size (as shown in figure 8a), the mean filament spacing does not impose any constraints on the eddy size that can penetrate the canopy.In this situation, the mean location of the virtual origin, y vo , is determined by the average distance between the eddy core and the impermeable bottom wall (Monti et al. 2020).Conversely, in a dense canopy regime, when S/h < 1, only eddies with sizes smaller than S can fit within the void regions of the canopy (see figure 8c).If an external eddy exceeds S, its interaction with the filaments leads to a local deformation, resulting in a penetration length smaller than S. In this situation, the virtual origin position, seen by the outer flow, moves closer to the canopy tip at a distance of approximately S from it.The mean location of y vo would scale directly with S (see figure 9), which supports the hypothesis proposed by MacDonald et al. (2018) for turbulent flows over spanwise aligned bars.The transitional state between the dense and sparse asymptotic conditions occurs when the canopy aspect ratio approaches unity ( S/h ≈ 1).In this regime, the void regions can be approximated as square cavities accommodating eddies with a size of ≈ S, resembling turbulent flows over d-type rough surfaces (Perry, Schofield & Joubert 1969;Leonardi, Orlandi & Antonia 2007;MacDonald et al. 2016).The mean location of the virtual origin in this case is positioned at the midpoint of the canopy height (as visible from figure 9), supporting the proposed transition criterion.
Figure 10 provides a qualitative assessment of the transition criterion, showing the contours of streamwise and wall-normal vorticity fluctuations on a y-z cross-section of the channel for different λ values.In the sparser regime, the outer flow fully penetrates the canopy layer, while in the denser regime, the penetration depth of the outer vortical structures is influenced by S. The localised fragmentation of large-scale vortical structures at the canopy interface in the densest cases aligns with the model predictions.By establishing the criterion for the onset of the dense regime ( S − d ≈ h), we can provide a precise estimate of the mean filament spacing and the corresponding solidity value (λ) that leads to the canopy flow transition.Considering the specific geometric parameters of the canopy used in this study (h/H = 0.1 and d/h = 0.24), we obtain S/H = 0.124,  corresponding to λ = (hd)/ S 2 = (H/ S) 2 (d/h)(h/H) ≈ 0.15 or S/h = 1 + d/h ≈ 1.24.The transitional λ value of 0.15 closely aligns with previous findings by Monti et al. (2020).It is also noticed that a similar threshold value of λ ≈ 0.15 is commonly accepted for distinguishing sparse from dense regimes in k-type rough surfaces (Schlichting 1937;MacDonald et al. 2016).

Statistical characterisation of the flow
We examine the mean velocity profiles by separately considering the region within the filamentous layer from the bed to the virtual origin and the outer region above y vo .Figure 11(a) presents the mean velocity distribution in a semi-logarithmic scale for all the solidity values.The profiles are normalised using two friction velocities: u τ,in = √ τ w /ρ at the impermeable wall (y = 0) and u τ,out = √ τ (y vo )/ρ at the virtual origin (y = y vo ), where τ w is the wall shear stress at the bed and τ (y vo ) is the total stress at the virtual origin.The figure shows that regardless of the solidity value, the mean velocity profiles collapse within the first few wall units ỹ+ (see figure caption for its definition), indicating the dominance of wall friction over the drag from filaments, even for denser canopies.(Hama 1954) plotted as a function of λ for all considered canopies.The solid line represents a polynomial fit for cases obtained varying S (current investigation), while the dashed line represents a polynomial best fit for the cases considered by Monti et al. (2020).See table 1 for the symbols' meaning.
As we move away from the wall, the modulation imposed by the filaments and the effects of the external flow start to play a role in the mean velocity distribution.
Outside the canopy, the mean velocity profile is scaled with the outer friction velocity u τ,out and the corresponding viscous length scale δ ν,out = u τ,out /ν.It follows a logarithmic distribution, and the presence of the canopy causes a shift in the logarithmic region (see inset in figure 11a).In this case, the logarithmic law for turbulent boundary layers over generic textured surfaces can be expressed as The first three terms in (3.4) represent the canonical log-law for a smooth wall, while the last term U + corresponds to the roughness function (Hama 1954;Jiménez 2004).The latter accounts for the momentum deficit or surplus due to the presence of textures on the wall and determines the wall offset in the logarithmic region.Examples of U + variations can be found in studies on permeable substrates (Rosti, Brandt & Pinelli 2018; Gómez-de Segura & García-Mayoral 2019) and rough wall boundary layers (Bechert, Bruse & Hage 2000;Jiménez 2004).
Figure 11(b) shows the roughness function as a function of the solidity λ for all the canopies considered.All the canopy configurations considered in this work exhibit a drag increment, with the magnitude of the roughness function depending on λ.A polynomial fit is provided for the cases where S is varied in the current investigation, while a dashed line represents the best fit for the cases considered by Monti et al. (2020).The symbols used in the figure are described in table 1.
Next, figure 12 and k + eff for the canopies considered by Monti et al. (2020) show a monotonous increase with λ, with a possible saturation of S + for very dense cases.However, when increasing the filament packing, a decrease in both S + and k + eff is observed for increasing values of λ.The effective canopy height and mean spacing for the densest case (i.e.case HD) are found to be S + ≈ 50 and k + eff ≈ 20, respectively, which correspond to the conditions on the protrusion height for the breakdown of the drag reduction regime in walls covered with riblets (Bechert et al. 2000;Garcia-Mayoral & Jimenez 2011).A decreasing effective height and drag in the dense regime (i.e. when λ > 0.15) can be related to the departure from a fully rough regime and the onset of a transitionally rough regime as previously reported by Jiménez (2004) and by Thakkar et al. (2018).The transitional regime is characterised by a complex dynamical interplay between the drag offered by the roughness elements and the viscous drag (Flack 2018).
The equivalent sand roughness k s is another parameter that can be used to characterise the effect of non-smooth walls (Schlichting 1937;Bradshaw 2000;Jiménez 2004).Here, k s was introduced as a notional measure of the height of closely packed sand grains that would produce the same frictional drag in the fully rough regime (when k + s ≥ 70) (Nikuradse et al. 1950).Based on the collapse of the roughness function in the fully rough asymptotic condition (Flack & Schultz 2014), k s can be defined via a modified law of the wall (Jiménez 2004) as The ratio of k s /k (where k is the roughness element height) and the roughness solidity are strongly related as shown in figure 13(a), where the non-dimensional form of the equivalent sand roughness (k s /k) is plotted as a function of the effective solidity λ eff = dk/ S 2 for all the canopy configurations considered here and by Monti et al. (2020).In the graphs, the value of k s /k has been rescaled with a filament drag coefficient C d (Jiménez 2004) computed at the mid-location of the filaments (the drag at any Lagrangian node along each filament can be readily obtained within the framework of the immersed boundary method we use (Pinelli et al. 2010) as the force per unit volume necessary to establish the zero-velocity condition).The figure indicates that all the considered canopies belong to the sparse k-type regime, with the canopy-dense cases (HD, DE, MDM and MDE) in a transitional condition.Figure 13(b) shows the relationship between k s / k and λ (i.e. the h based canopy solidity).In the sparse k-type regime, k s / k increases linearly with  (Jiménez 2004) computed at the mid-location of the filaments, where the localised flow is unaffected by the fixed and free ends.The dashed lines are from Jiménez (2004).See table 1 for symbols.Other symbols represents non-circular roughness elements: * , spanwise fences (Schlichting 1937); 'x', spanwise fences (Webb, Eckert & Goldstein 1971); , spanwise cylinders (Tani 1987).λ, reaching a maximum at the transition threshold (λ = 0.15).In denser regimes, the value of k s / k decreases with increasing λ due to the sheltering effect of the filaments.While the concept of equivalent sand roughness can be helpful in establishing a connection between the outer behaviour of canopy flow and traditional experiments conducted on rough walls, it is evident from our findings that solely examining the external canopy flow is inadequate for fully characterising the overall impact of the canopy on the flow development above the canopy tip.In fact, the flow behaviour above the canopy exhibits similarities to various roughness morphologies, which are related to the solidity parameter.
The diagonal Reynolds stresses are analysed and presented in figure 14, where they are scaled with the outer friction velocity.The profiles collapse in the region outside the canopy, while the inner layer exhibits a substantial influence of the in-plane solidity in all three components (streamwise, wall-normal and spanwise).The outer region profiles collapse similarly to smooth wall turbulent channel flows, confirming outer layer similarity.The spanwise component distributions show an increase in the magnitude of the internal peak moving towards the canopy bed as the canopy solidity increases.The streamwise component variations exhibit a non-monotonic behaviour, with a decrease in magnitude for sparser regimes and an abrupt increase for denser regimes.The increase is accompanied by a movement of the peak towards the canopy bed, indicating the presence of another mechanism associated with the onset of the dense regime.
The diagonal Reynolds stress profiles are further analysed using a local velocity scale based on the total stress: This is the same scaling proposed by Sharma & García-Mayoral (2020a) in the context of canopy flows, and by Jiménez (2013) and Tuerke & Jiménez (2013) when dealing with manipulated wall-bounded turbulent flows.The profiles scaled with this local friction velocity (presented in figure 15) show a collapse across all cases, resembling smooth-wall-like behaviour, particularly in the canopy region.They exhibit an almost bimodal distribution, with two peaks located in the inner and outer regions.The magnitude of the inner peak increases with λ for all three components.A detailed examination of the wall-normal component reveals that while the peak magnitude increases with λ in the inner layer, it remains unaffected by the stem number density.The spanwise component profiles show an increase in the magnitude of the internal peak with increasing λ.The streamwise component variations exhibit an interesting behaviour, with a decrease in magnitude for sparser regimes and an increase for denser regimes, indicating the influence of the canopy density.

Structure of the flows
We now turn our attention to the inception and arrangement of coherent structures within both the intra-canopy and outer regions, with a particular focus on the role of the canopy's in-plane solidity.We begin by examining the time-averaged flow field using a spatial ensemble averaging technique.Subsequently, we delve into the fluctuating flow field through triple decomposition and investigate the energy content of velocity fluctuations using spectral decomposition and 3-D autocorrelation functions.
The time averaged flow fields are obtained following a three-stages procedure (Monti et al. 2020).First, we compute a time-averaged field encompassing the entire canopy.This mean field is not particularly meaningful as it contains local variations resulting from the random distribution of filaments across the tiles.To mitigate this effect, we proceed to the second stage, where we introduce a virtual cuboid with a base measuring 2 S × 2 S and a height equal to h.In this step, we adjust the velocity fields over each tile volume by shifting them to align their respective stems with the centre of the cuboid base.All the resulting translated fields, occupying the cuboid, are then ensemble-averaged to generate an intermediate mean field.Finally, in the third stage, we obtain the double-averaged field over a S × S tile, akin to what would correspond to a uniform filament distribution, by averaging across the x and z directions of the cuboid.This is achieved by leveraging the averaged velocity field's periodicity and parity.
Figure 16 presents the mean flow fields obtained using this procedure, displaying the mean streamlines and iso-contours of the wall-normal velocity.For visual clarity, the filaments located at (x/ S, z/ S) = (−1, 0) and (0, −1) in the HD configuration have been removed.In the sparser regime (figure 16a,b), the plane cut at the canopy tip reveals a mean ejection and sweep on the frontal and rear sides of the filaments.However, deeper within the canopy, we observe an opposite behaviour with the mean ejection and sweep occurring on the rear and frontal sides of the filaments at y/h ≈ 0.8, indicating that the outer flow penetrates the canopy by sweeping along the rear side of each filament.Closer to the wall, the mean flow pattern becomes independent of the wall-normal coordinate.Furthermore, it is observed that the footprint of the mean sweeps observed near the canopy tip progressively diminishes as the canopy density increases (as seen in figure 16b-d).The sparser configurations, especially those in the mildly dense regime (i.e.MD and DE), exhibit very weak sweeps on the rear sides of the filaments, suggesting that the mean flow field is nearly unaffected by changes in the wall-normal coordinate.
An interesting observation arises from the mean flow analysis of the HD configuration (see figure 16c).It reveals the presence of a vortex between adjacent filaments, suggesting a decoupling phenomenon within denser canopy regions.This decoupling inhibits direct kinematic communication between the inner and outer flows, resulting in the emergence of a vortex with a size determined by filament spacing and local Reynolds number.This phenomenon bears resemblance to observations made in previous studies on flows over k-type rough surfaces (Perry et al. 1969;Jiménez 2004;Leonardi et al. 2007), further emphasising the similarity between the present HD configuration and these rough surface conditions.
In the mildly dense regime (MD and DE configurations), the mean flows resemble those of the fully dense configuration (MDE) studied by Monti et al. (2020).Partial decoupling is observed between the internal and external flows, limiting the formation of a stable vortex.The emergence of a mean vortex in the HD configuration can be conceptualised as a cavity-like flow, following principles proposed by Prasad & Koseff (1989) and Hou et al. (1995).This similarity becomes evident by looking at figure 17(a) that reports the rescaled, averaged streamwise and wall-normal velocities using the internal friction velocity (u τ,in = (τ w /ρ) 1/2 ) for the densest case within the void between successive filaments.The velocity distribution of the streamwise velocity components shown in figure 17(a) also highlights the presence of a sheltering effect from the upstream filament.This effect (Poggi et al. 2004) is confirmed by the mean streamlines on the wall-parallel plane near the canopy tip (figure 16).In the sparser configurations, the streamlines show deflections around the filaments rejoining downstream.Differently, in denser canopies, the mean streamlines on the rear side of filaments do not reconnect, indicating a stronger sheltering effect.This is particularly evident in the HD case.Increased sheltering leads to reduced U + (figure 11b), as observed in studies on rough surfaces with densely packed elements (Jiménez 2004).
We now shift our attention to investigating a potential instability within the canopy, as suggested by Monti et al. (2020).To explore this phenomenon, we examine figure 18, which presents instantaneous snapshots of velocity fluctuations in wall-parallel planes at locations corresponding to internal inflection points.In sparser configurations (SP and TR), the velocity fluctuations exhibit coherent regions resembling external streamwise-oriented vortices, along with smaller-scale fluctuations induced by fluid meandering around the filaments.However, in denser canopies, more pronounced modulation in the internal flow field becomes apparent.Panels ( j) and (m) highlight strong spanwise-oriented structures resembling streamwise-moving fronts.Monti et al. (2020) proposed that these spanwise coherent motions, reminiscent of Kelvin-Helmholtz-like instabilities, arise from high wall-normal permeability that promotes the directional penetration of external structures.Within the canopy, these penetrations generate streamwise and spanwise modulations due to bed impermeability and solenoidal conditions.
The variation and organisation of coherent motions within the canopy layer can be quantified through higher-order statistics of the full fluctuating field, specifically, the two-dimensional two-point auto-correlation function in the streamwise ( x) and spanwise ( z) directions, yields the structure of the velocity fields associated with various regimes.Figure 19 illustrates the positive (0.1) and negative (−0.1) correlation coefficients, showing x/H  alterations in the red iso-surface related to the wall-normal velocity component correlation, except for the HD configuration, which exhibits a spiky structure.Conversely, the negatively correlated iso-surface (blue region) shrinks for larger λ values until disappearing in the densest case.
The streamwise correlation function reveals that an increase in the in-plane solidity modifies the coherence length and direction, while the SP case shows long-scale correlation above and across the canopy.This behaviour links to the quasi-streamwise vortices that penetrate the canopy through sweep events.
As the canopy density increases, the penetration of external vortices into the canopy is governed by S. The HD configuration shows the shortest streamwise coherence length at the canopy top, where the vortices lose their coherence due to deformation.Furthermore, streamwise coherence decreases with proximity to the bottom wall, while spanwise coherence becomes more prominent.
The fluctuating flow field structure can be further explored by considering the spectral energy content of the fluctuating velocity components.The one-dimensional premultiplied spectra, rescaled with the local friction velocity, are shown in figures 20 and 21.They feature two peaks corresponding to the wavelengths λ x /H 1 and λ z /H 1, which are connected to the coherent structures populating the external flow.
The internal spectral peak (y/H ≤ 0.1) points towards two scenarios arising in different canopy flow regimes.In the sparser regimes, the inner peaks do not correspond to the canopy geometrical parameters' wavelengths.For larger λ, a peak emerges corresponding to the filament's diameter.This peak is most intense for the HD configuration, highlighting the wake's significance behind the filament in the dense regime.
In denser configurations, the internal peak of the spectral energy content splits into two maxima.The leftmost peak corresponds to the average filament spacing (λ x /H = λ z /H = S), while the rightmost mirrors the outer peak.Finally, Monti et al. (2020) associated the peaks of u and w with the energetic sweeps of the outer layer vortices into the canopy and subsequent deflection upon reaching the bottom wall, affected by the wall's impermeability and solenoidal constraints.The HD case exhibits a different behaviour due to the increased filament clustering affecting the wall-normal permeability, thus reducing the outer vortices' role in determining the intra-canopy region's dynamics.
The internal structure of the intra-canopy layer is elucidated by considering the spectral decomposition of the fluctuating flow field.Figure 22 displays the two-dimensional premultiplied spectra in log coordinates, differentiated by colours.The red, blue and yellow regions denote the inner region, outer layer and overlap zone, respectively.Red and green lines illustrate wavelengths linked to the canopy height and mean filament spacing.
In sparser conditions (panels a-c), iso-surfaces span the canopy, exhibiting a wide range of wavelengths from λ x /H = λ z /H = 1 to those tied to canopy geometry.Distinctive peaks within the spectra can be attributed to large-scale vortices infiltrating the canopy layer and smaller contributions from meandering filament motion.
With increased canopy density, the spectra's iso-surfaces reshape more intricately.As the solidity λ rises, there is a more defined scale separation and energy content changes based on the wall-normal distance from the canopy bed.An elongated peak appears within a wavelength window of λ x /H, λ z /H ∈ [h, S] across all components of the premultiplied spectra, associating with velocity fluctuations due to filament interaction.
Interestingly, the v component's spectral footprint displays two elongated peaks.Peaks extend further with increasing λ values.A mildly elongated peak is noticed in the u component in the TR case (panel d), possibly due to transitional effects present at this λ value.
In the HD case, a pronounced scale separation occurs in the intra-canopy region, visible in panels (m,n,o).The u component exhibits three isolated peaks at varying wavelengths, while the v and w components display peaks that correlate to largeand small-scale contributions.This suggests a correlation between the three velocity fluctuation components within the canopy.Peaks of the v component progressively decrease with an increasing λ. representing the premultiplied spectra κ x Φ u u /u 2 τ,l , κ x Φ v v /u 2 τ,l and κ x Φ w w /u 2 τ,l , respectively.Solidity values increase moving from top to bottom.The yellow, red and green vertical lines correspond to the wavelengths matching the size of the diameter, the filament height and the spacing respectively.The cyan, yellow and red horizontal dashed lines represent the wall-normal location of the virtual origin, the location of the internal inflection point and the canopy tip.One-dimensional premultiplied spectra do not clearly identify coherent internal motions.Large-scale contributions are confined to the canopy's upper region in the HD case, while small-scale contributions reside in the wavelength window λ Finally, our analysis returns to the external boundary layer above the canopy.Turbulence in this region is known to provoke large-scale spanwise coherent structures, identified either through the two-dimensional spectral content of wall-normal velocity fluctuations above the canopy tip, or the outer peak of the premultiplied cospectra of u and v The Kelvin-Helmholtz-type instability, analogous to a plane mixing layer, arises due to an inflected mean velocity profile.Many have used this analogy to study coherent structures at the canopy tip (Raupach et al. 1996;Finnigan 2000;Nepf 2012).Authors hypothesised that the spanwise roller wavelength (Λ x ) in a dense canopy, scaled with the shear length (L s ), falls within 7 < Λ x /L s < 10.Raupach et al. (1996) proposed a precise estimate of Λ x 8.1L s .from the tip, which does not account for shorter canopies.To include sparser canopies, we have adjusted the calculation to U o = U h − L s (dU/dy) y=h using the mean velocity at the virtual origin (U h being the velocity at the canopy tip, i.e. at y = h).
Our findings, presented in figure 25(b), indicate that shear length scale's validity as vorticity thickness is constrained to canopies with large velocity differences, specifically those exceeding a threshold magnitude of U s /U b ≈ 0.26.

Conclusions
We have conducted a series of detailed numerical simulations of turbulent open channel flows developing over submerged canopies composed of rigid, slender, circular cylindrical stems mounted perpendicular to an impermeable wall.Under these conditions, the hydrodynamic regime is principally determined by two geometrical parameters: the height of the stems (h) and their mean spacing ( S).
For slender filaments (i.e.filaments with a diameter d such that d/h 1), the non-dimensional combination of the canopy geometrical parameters λ = dh/ S 2 (referred to as the solidity) is generally deemed a sufficient indicator, enabling the prediction of the flow regime within and above the canopy.A previous, companion study by Monti et al. (2020) examined the impact of the canopy height on the onset of various regimes without altering the average spacing S between the stems.In contrast, this study has adopted the complementary approach of changing the value of λ by keeping the height of the stems constant while adjusting the number density of the canopy (i.e.modifying the value of S).This work, combined with the findings of Monti et al. (2020), has led to the following primary results: (i) the development of a conceptual model encapsulating the core elements of the transition from a sparse to a dense canopy regime (Nepf 2012); (ii) an understanding of the role played by λ, h and S in determining the nature of the intra-canopy flow; (iii) the identification of the relevant length scales in the velocity layers characterising the canopy flow; (iv) the provision of a detailed description of the turbulent flow fields within and above the canopy under various regimes, using a triple decomposition technique.
Broadly speaking, this study has found that sparse and mildly dense configurations (i.e.λ 0.26) share several similarities with canopies of the same solidity obtained by varying the height of the filament (Monti et al. 2020).However, under denser conditions, the overall intra-canopy flow experiences a slowdown due to a decrease in the streamwise and spanwise permeabilities.
When comparing the highest solidity cases of the current study with those reported by Monti et al. (2020) (which share the same magnitude of λ = 0.56), it is evident that the flow dynamics in the canopy region differ.In the case of tall canopies considered by Monti et al. (2020), the transfer of momentum between the canopy's outer and inner layers occurs in the form of wall-normal sweeps, driven by the external flow.When these jet-like structures approach the wall, they are deflected in the streamwise and spanwise directions, re-energising the otherwise stagnant flow in the vicinity of the wall.
In the densest case, obtained by clustering the filaments, the redistribution of momentum in the near-wall region is partially inhibited by the decreased permeability parallel to the wall, which results from increased filament packing.Under this condition, the intensity of the wall-normal momentum directed towards the wall is weaker than that observed by Monti et al. (2020).This momentum deficit also prevents the formation of the kink in the mean velocity profile observed by Monti et al. (2020) in their densest (i.e.highest λ) case.
The configurations considered in this study have confirmed that the solidity value λ defines the sparse and dense, asymptotic regimes.
The transition from the sparse to the dense regime is governed by the canopy aspect ratio, defined as the ratio between the mean filament spacing and the filament height.In the transitional condition (i.e.S/h = 1), our simulations have demonstrated that the mean location of the internal inflection point and the virtual origin coincide at a single point (y/h ≈ 0.5).We have shown that this condition corresponds to λ = 0.15 ≈ 0.2, which is very close to the commonly accepted value that separates sparse and dense regimes in the context of canopy flows (Nepf 2012;Monti et al. 2020) and flows over rough surfaces (Schlichting 1937;MacDonald et al. 2016).
At the onset of the dense condition, the parameter S begins to play a leading role in controlling the size of the eddies that can cross the canopy/external flow interface by setting the effective filter width for the external structures.By determining the size of the eddies permitted to penetrate the canopy, the mean filament spacing S essentially sets the location of the virtual origin observed by the outer flow.
We also found that the drag exerted by the canopy on the outer flow (i.e.U + ) starts to decrease when more densely packed canopies are considered.This behaviour can be attributed to the mutual sheltering effect of the downstream filaments, which becomes significant in the dense regime.
This research has also facilitated the identification of the length scales that best characterise the flow layers developing inside the canopy.Making the thickness of the internal layers non-dimensional with the appropriate length scales enables the unravelling of universal features of the flow inside the canopy.In particular, we have established that the internal layer scales well with S, suggesting that the nature of the inner flow is solely governed by the mean filament spacing, regardless of how λ was varied.
In the S → 0 limit, the internal layer thickness must collapse to zero since the wall-normal sweeps are entirely inhibited for very low stem spacing.Conversely, when the solidity is increased by considering taller canopies, the inner layer thickness becomes asymptotically independent of h.
Regarding the external canopy layer, it was found that the penetration depth of the outer flow is mainly set by the canopy height h, once again independent of the overall value assigned to λ.However, since the mean filament spacing S controls the size of the eddies allowed to enter the canopy, the level of penetration decreases to zero as S → 0 to satisfy the solid wall condition.Conversely, when the canopy height is increased towards the dense regime, the level of penetration saturates to S.
To identify the coherent structures populating the canopy and the external regions (figure 26 provides a qualitative overview of the structures within the flow), we used spectral analysis and the velocity autocorrelation functions.As highlighted by many other authors, our work also found that the inflection point located at the tip of the canopy (arising as a result of the abrupt drag discontinuity) triggers an instability that generates localised Kelvin-Helmholtz-like rollers enhancing the aforementioned wall-normal momentum transfer into the canopy.
The occurrence of these spanwise coherent structures has been previously reported by many authors in the context of different shear flows over textured surfaces: in dense canopy flows (Monti et al. 2020;Sharma & García-Mayoral 2020b), above permeable substrates (Rosti et al. 2018;Gómez-de Segura & García-Mayoral 2019) and ribletted surfaces (Garcia-Mayoral & Jimenez 2011).the momentum directed towards the wall is redistributed; a fraction of it is repelled from the wall, giving rise to regions of high Reynolds stresses (i.e.u v peaks).
The combination of these mechanisms, primarily driven by the transfer of momentum towards the wall, essentially sets the location of the internal inflection of the mean velocity profile (Monti et al. 2022).The outer structures embedded in the backdrop of the outer logarithmic region, capable of penetrating the canopy in the transitional and mildly dense regimes, also leave their footprints within the canopy layer.However, the memory of their presence diminishes as the solidity is increased, eventually vanishing in the fully dense regime.In this condition, the inner coherent motions are almost completely isolated from the outer flow, and a decoupling effect takes place, producing a large-scale separation between the inner and the outer layers in the canopy.

Figure 2 .
Figure 2. Computational box of the current study.
Figure 3(a) presents the dimensionless mean velocity distributions normalised by the open channel height (H) and the bulk velocity (U b = (1/H) H 0 U dy, where • represents the triple averaging operator encompassing time and homogeneous directions).

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

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

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

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

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

Figure 11
Figure 11.(a) Mean velocity profiles normalised with the two friction velocities u τ,in = √ τ w /ρ and u τ,out = √ τ (y vo )/ρ (see text).The abscissa ỹ+ represents the wall-normal coordinates in wall units determined using either the inner (left set of profiles for y + in < 5) or the outer friction velocity (right set of profiles for y + out > 50) and considering either the wall or the virtual origin as the zero of the wall-normal axis: ỹ+ = yu τ,in /ν and ỹ+ = (y − y vo )u τ,out /ν.The solid line without symbols refers to the profile of a canonical smooth wall channel flow at Re τ = 950 (Hoyas & Jiménez 2008).The sub-plot in panel (a) is a magnification of the profiles in the logarithmic region.(b) Roughness function(Hama 1954) plotted as a function of λ for all considered canopies.The solid line represents a polynomial fit for cases obtained varying S (current investigation), while the dashed line represents a polynomial best fit for the cases considered byMonti et al. (2020).See table 1 for the symbols' meaning.

Figure 12
Figure 12.(a) Mean filament spacing and (b) effective canopy height expressed in outer wall units (u τ,out ), as a function of λ for all considered canopies.See table 1 for symbols.

Figure 13 .
Figure 13.Non-dimensional form of equivalent sand roughness (k s /k) as a function of the (a) effective solidity (λ eff ) and (b) solidity (λ).The value of k s /k has been rescaled with a filament drag coefficient C d(Jiménez 2004) computed at the mid-location of the filaments, where the localised flow is unaffected by the fixed and free ends.The dashed lines are fromJiménez (2004).See table 1 for symbols.Other symbols represents non-circular roughness elements: * , spanwise fences(Schlichting 1937); 'x', spanwise fences(Webb, Eckert & Goldstein 1971); , spanwise cylinders(Tani 1987).

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

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

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

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

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

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

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

Figure 21 .
Figure 21.Premultiplied spectra of the fluctuating velocity components κ z Φ u u /u 2 τ,l , κ z Φ v v /u 2 τ,l and κ z Φ w w /u 2 τ,l represented as a function of the spanwise wavelength (λ z /H) and wall-normal distance (y/H).For the panels' ordering and meaning of the coloured lines, refer to the caption of figure 20.

Figure 22 .
Figure22.Iso-surfaces of two-dimensional premultiplied spectra of the fluctuating velocity components represented as a function of the wall-parallel wavelengths (λ x /H and λ z /H) and the wall-normal distance (y/H).The first column (a,d,g, j,m), the second column (b,e,h,k,n) and third column (c, f,i,l,o) represent the κ x κ z Φ u u /u 2 τ,l , κ x κ z Φ v v /u 2 τ,l and κ x κ z Φ w w /u 2 τ,l premultiplied spectra, respectively.The rows from top to bottom correspond to decreasing in-plane solidity.See the text for the meaning of the lines and colours.Note that the isosurfaces are open-ended in the wall-normal direction because the region above the canopy is not reported here.

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

Table 1 .
Monti et al. (2020)rs considered in both the current investigation and in the one byMonti et al. (2020)(cases identified with an * in the table): λ is the solidity defined in (1.1); S/h is the canopy aspect ratio; n i × n k is the filaments cross-plane count, where n i and n k are the numbers of filament rows in the streamwise and spanwise directions, respectively; Re τ,o is the friction Reynolds number based on the friction velocity (u τ,o