Turbulent structure inside and above shallow to deep canopies

Multi-plane PIV measurements were performed in an open-channel flume filled with elongated prisms of height k and width to investigate the effect of the deepening of the canopy on the flow structure. Velocity measurements were performed both inside the canopy and above it. Analysis of the spatial convergence for the double-averaged quantities shows that for canopy flow investigations (z < k), at least 5 measurement planes are required to obtain a relative spatial convergence error below 3% for the dispersive shear stress, the quantity the most sensible to spatial sampling. With only three measurement planes, the spatial convergence is below 1% only in the flow region above the canopy (z > k). Three canopy aspect ratios, k/ = [1, 3, 6] were investigated for a fixed modified-submergence ratio β = (h − k)/ = 3 where h is the water depth. As the canopy deepens, the hydraulic roughness decreases and the velocity near the bottom of the canopy becomes gradually constant, as expected for deep canopies. We show how the highly converged (both in space and time) profiles of double-averaged longitudinal velocity and total shear stress can be used to calculate the vertical distribution of drag in the canopy. With this methodology, values of the drag coefficient CD(z) can be calculated, and are found to be always close to unity, even in the upper part of the canopy.


Introduction
Turbulent boundary layer flows in and above urban or plant canopies have been studied extensively, for roughness elements geometry ranging from cubes [4] to elongated cylinders [5].In 2007, Nikora et al. [6] showed the relevance of the double-averaging methodology to describe these flows in the roughness sublayer where time-averaged fields remain fully three-dimensional.Above the canopy (for z > k where k is the height of the roughness elements), the mean velocity field is well known for high relative submergence flow conditions (i.e., high values of the relative submergence ratio h/k where h is the water depth, typically h/k > 15) and exhibits the classical logarithmic law behaviour.This logarithmic law was shown by Keulegan in 1938 [3] for rough or wavy surfaces.The existence of the logarithmic law for very low relative submergence ratios h/k, approaching unity, was recently shown to prevail down to h/k = 3 in [2] and even to h/k = 1.5 in [9] for a square arrays of cubes.
By analysing velocity measurements performed above (z > k) and inside the canopy (z < k), Nepf [5] distinguishes two kinds of canopies: sparse canopies for which the resistance of the bed of the canopy is as important as the drag of canopy elements and the velocity profiles remain logarithmic through the canopy, and dense canopies, for which the bed resistance becomes small in comparison with drag exerted by the canopy elements, leading to the formation of an intense shear zone and turbulence source near the top of the canopy.
To observe the effect of the canopy depth (i.e. the frontal aspect ratio of the roughness elements composing the canopy) on the turbulent structure, within and above the canopy, three configurations were investigated and are discussed in this paper.They correspond to the same modified-relative-submergence ratio defined here as β = (h − k)/ (where is the width of the roughnesses).The first configuration is the same cubic rough-bed configuration as the one named S1 in Florens et al. [1] and Eiff et al. [2], with an aspect ratio k/ = 1 modelling urban zones.The two others configurations make use of elongated prisms, representing vegetation (or high buildings), with aspect ratios k/ = 3 and 6.
To investigate the velocity field inside (0 < z < k) and above (z > k) the canopy, PIV measurements (Particle Image Velocimetry) were performed in several planes of a typical roughness pattern.Double-averaged quantities could then be inferred from this measurements and used to estimate the drag distribution in the canopy.The experimental set-up is presented in section 2. The accuracy of the double-averaging methodology for a finite set of measurement planes is discussed in section 3. Analysis of the turbulent structure and drag distribution is done in section 4.

Experimental set-up
The experiments were performed at the Institut de Mécanique des Fluides de Toulouse (IMFT) in a 26-m-long, 1.10-m-wide and 0.50-m-deep open-channel flume with a slope of 0.3 %.To obtain uniform flow conditions, the water level was adjusted with a downstream weir.The water discharge ranges from 1 L.s −1 up to 150 L.s −1 , with an accuracy of ± 0.2 %.Within the entry section, the low passes through a honeycomb and a series of mesh grids before converging at the main channel entry to help establish inlet low homogeneity and to reduce both water surface waves and the background turbulence intensity.
The bottom of the open-channel was covered with prismatic roughness elements (k × × , where k and are respectively the roughness element height and width.Different rough beds were investigated for canopies ranging from the urban type (k/ = 1) to the vegetation type (k/ = 6).The first rough bed, with an aspect ratio k/ = 1, is a square distribution of 4 cm cubes with a planar density λ p of 0.2 (λ p = 2 /L 2 where L is the distance between the rough elements).The second rough bed was built with elongated prisms (12 cm × 2 cm × 2 cm) with an aspect ratio k/ = 6, conserving the same planar density.An intermediate canopy configuration was also investigated with smaller prisms and an aspect ratio k/ = 3 (6 cm × 2 cm ×2 cm), again with the same planar density.The first rough-bed (k/ = 1) extends longitudinally over 10 m and the two others extend over 5 m, in order to obtain fully developed boundary layers in the measurement section, located at x M = 19.20 m (where x M is the center of the measurement area).For all these bed geometries the modifiedrelative-submergence ratio β was kept constant with β = 3.1.In order to investigate only the effect of the canopy depth.Table 1 summarizes the experimental parameters for the three configurations.Table 1.Flow regimes investigated.h is the roughness element height, the roughness element width, k/ the roughness aspect ratio, λ p the planar density, λ f the frontal density (λ f = h /L 2 ), Q the water discharge, D the water depth, U bulk the bulk velocity, D/h the relative submergence, β = (h − k)/ the modified-relative-submergence.To allow optical access for Particle Image Velocimetry (PIV), BK7 glass roughness elements were used in the measurement section.Using a 2 × 200 mJ pulsed Nd-Yag connected to a parallel-laser-sheet generator, negligible shadows were produced when the laser sheet passed through the flume bottom and glass roughness elements.Instantaneous twocomponent velocity fields were measured in at least five parallel vertical planes uniformly distributed along one periodic pattern with the locations: y = /2, 3/4 , , 1/4(L − ) and 1/2(L − ) (see figure 1), where y is the transversal direction to the flume.The laser sheet thickness was 2 mm.In all planes, the images were recorded through the glass side-wall of the flume with a high resolution camera (2k×2k pixels) and an acquisition frequency of 3 Hz.A telecentric lens was used to prevent any parallax effects since the optical path includes the transparent glass prisms in water.All PIV measurements were performed near the center of the flume with 10 000 decorrelated velocity fields and a spatial resolution of 0.6 mm.This set-up allowed measurements along the whole flow depth, from z = 0 to z = h (position of the free surface).All PIV post-processing was done with the PIV software CPIV-IMFT developed at IMFT (a FFT cross-correlation method with peak-locking reduction schemes and parallelization).

Double-averaging methodology
The roughness-layer is defined as the region of the flow with spatial inhomogeneity for the mean (time-averaged) quantities, and extends above the canopy top.In this flow region, it is relevant to adopt a double-averaging methodology as introduced and discussed by Nikora et al. [6] for free-surface flows.It consists in a decomposition quite similar to the Reynolds decomposition, where an instantaneous quantity θ can be decomposed into : where θ i and θ are respectively the temporal and the spatial fluctuations.Using the double-average methodology in the Navier-Stokes and continuity equations and assuming a two-dimensional, steady and uniform flow over a rough-bed, the longitudinal momentum and continuity equations read : where g is the gravitational acceleration, γ the slope of the channel, F x is the total drag force (including both viscous and wake drag forces) and Φ is the porosity defined as: Φ = Φ(z) = S f /S t , where S f and S t are, respectively, the fluid surface and the total surface (fluid and solid) at level z.The total shear stress τ i j is defined as the sum of the dispersive stress, the spatially-averaged Reynolds shear stress and the spatially-averaged viscous stress: Above the canopy, without obstacles, there are not drag forces and Φ = 1 so that equation 2 can be integrated from z = h (where the total shear stress is equal to zero) to z, in order to obtain a total shear stress that reads τ xz (z) = ρg(k−z) sin γ.Within the canopy, it is possible to keep on integrating equation ( 2) from z = k, taking into account the drag forces and assuming continuity of φτ xz at z = k.There, the total shear stress reads : This last equation, when applied at z = 0, shows that the total drag on the flow (the sum of the bottom drag τ xz (z = 0) and the drag along the canopy given by h 0 F x dz) is equal to g sin (γ) (h − k + Φk), as proposed by Pokrajac et al. [7].The profiles of τ xz and the forcing term g sin (γ), it is then possible to calculate the distribution of F x (z), as will be done in section 4.

Spatial convergence of the double-averaged flow quantities
In most existing experiments, only few vertical profiles have been used to estimate the doubleaverage quantities defined by Nikora et al. [6], resulting in a low level of spatial convergence (as pointed out by Florens et al. [1]).Here, to obtain spatially converged double-averaged quantities ( θ (z)), a discrete weighted average of each time-averaged field was used, following Florens et al. [1], and reading : where n is the number of vertical velocity profiles in a PIV slice and m the number of slices.δ i j and ∆x are, respectively, the width and the length of surface areas.
To investigate the influence of the number of PIV slices on the double-averaged quantities, we consider here that the estimated double-averaged quantities for the flow at k/ = 1 with nine different PIV measurement planes are spatially converged.Then, we subsample these measurements, using only three planes as in Florens et al. [1] and Eiff et al. [2], or five planes for the other configurations with higher values of k/ .Here, we define an absolute spatial convergence error for the quantity θ as θ = θ est − θ re f , where θ est is the estimated double-averaged quantity from three or five measurement planes and θ re f is the reference double-averaged quantity obtained with nine planes.
These absolute errors are normalized using the mainstream velocity u for velocity profiles or the friction velocity u * defined as u * = τ k /ρ (with τ k the total shear stress at the top of the canopy (z = k)) for second-order velocity tensor quantities.Figure 2 shows normalized spatial convergence errors for the double-averaged velocity components, Reynolds and dispersive shear stresses with five and three PIV slices.
Above the canopy (z/k > 1), u / u and w / u are very close to 0. This confirms that three planes are sufficient for the estimation of these double-averaged quantities above the canopy, as already shown by Florens et al. [1] by a comparison with measurements in horizontal planes.For second-order quantities, the convergence error is larger, even above the canopy, with errors as large as 4.5% for the dispersive shear-stress with only three planes.With lower sampling as is usual, considerably higher levels of error are expected for the dispersive stress, whose contribution to the total shear stress is important, which explains in part at least why in previous experimental studies, linear profiles of τ xz are scantly obtained down to the top of the rough bed for uniform flow conditions.
For the flow inside the canopy, with the three-plane estimates, there is an absolute error as large as 6.3% for the longitudinal double-averaged velocity estimate and the normalized error on the dispersive shear-stress is around 7.4%.With five planes, only a 3% normalized error remains for the dispersive shear-stress in the upper part of the canopy.Therefore, it was decided that five planes were required in order to adequately describe the flow within the canopy and all following double-averaged quantities were obtained with five planes.4 Turbulent structure

Mean flow and turbulence statistics
To observe the effect of the canopy depth we focused on the behaviour of the double-averaged mean and turbulence statistics for a constant modified-relative-submergence β = 3.1.Figure 3 shows the longitudinal velocity u and the total shear stress τ xz for the three different canopies (k/ = 1, 3 and 6).The almost perfectly linear behaviour of the total shear stress above the roughness elements (z > k) is in agreement with the assumption of uniform two dimensional and fully-developed flow conditions.For all flow regimes, the longitudinal doubleaveraged velocity above the canopy ((z − k)/ > 0) is well described by a logarithmic law and a wake law reading : where κ is the von Kármán constant (κ = 0.41), d is the displacement height, z 0 the roughness length, W(η) = (2Π/κ) sin 2 (η) the wake law with η = (z − d)/(h − d) and Π the Cole's wake parameter.
Deep inside the canopy with elongated prisms (with aspect ratio h/ = 6), the longitudinal velocity decreases towards a constant value.This trend is still observed for the intermediate canopy (with aspect ratio k/ = 3), but not at all for the shallow canopy with aspect ratio k/ = 1.For this last flow regime, an exponential decay is observed, as modelled by McDonald [4].The constant velocity found in the lower part of the deep canopy does not correspond to an equilibrium between the gravity forcing term and the wake behind the prisms.Indeed, in figure 3(right) there is a clear linear decrease of the total shear stress in the range −5.5 < (z − h)/ < −1.5 for aspect ratio k/ = 6.In equation (3), both the gravity forcing and the gradient of the total shear stress are constant and add up to yield a constant drag force F x .This constant drag, along with a constant drag coefficient assumption, yields the observed constant velocity.
All the log-law parameters have been calculated using the same methodology as in [2] and are given in table 2 which propose different normalizations.Moreover, the values of the loglaw parameters and log-law validity range for the shallow canopy with aspect ratio k/ = 1 are in agreement with values found in [2] for the same canopy and flow regime (α = h/k = 0.24 in Eiff et al. [2], corresponding to a modified-submergence aspect ratio β = 3.1).As k/ increases, the position of the log-law origin, given by d, moves deeper into the canopy, as shown by increasing values of (k − d)/ .Indeed, the displacement height d scales as k, with values around d/k = 0.8, almost constant for the three canopies investigated here (which have all the same value of λ p ).The hydraulic roughness k s , when scaled with k, exhibits first an increase up to k s /k = 0.91 for the intermediate canopy before decreasing again.The validity range of the log-law behaviour begins very close to the top of the roughness elements, and occupies almost 25% of the free flow region (k < z < h) for the intermediate and deep canopies.
Table 2. Logarithmic law parameters.u * , d, z 0 and Π are, respectively, the friction velocity, the displacement height, the roughness length and the Cole's parameter.z m and z M are the lower and upper bounds of the validity range for this log-law behaviour.

Total drag force and drag coefficient
Using equation ( 2), the total drag force within the canopy can be calculated as

Conclusion
The evolution of a canopy flow when the canopy depth increases with constant planar density (λ p = 0.2) and constant modified-submergence ratio (β = (h − k)/ =3 .1)was investigated with multi-plane PIV measurements.Using telecentric optical arrangements, the complete interstitial flow could be measured.It was shown that five planes are required to adequately estimate the double-averaged flow quantities, in particular the significant dispersive shearstress inside the canopy within 3% relative error.We also showed how these accurate doubleaveraged flow quantities can be used in uniform flow conditions to calculate the vertical distribution of the drag force inside the canopy.
For an elongated canopy (k/ = 6), deep inside the canopy, the flow behaves as a "pure" wake flow, without interaction with the top of the canopy and the canopy bed.The wake flow corresponds to the region where the longitudinal double-averaged velocity u is constant, as well as the drag force.This region still exists for the intermediate canopy with k/ = 3, but its vertical extent is reduced.For the shallow canopy with k/ = 1, this region disappears and the canopy-top mixing layer interacts directly with the canopy-bed boundary layer.Both the longitudinal double-averaged velocity and the total drag force profiles become more complex.
In future work, other free-surface flow conditions will be also investigated, corresponding to different values of the modified-submergence ratio aspect β, in order to study the combined effect of confinement by the free surface (decreasing β) and by the canopy bottom (decreasing k/ ).

Figure 1 .
Figure 1.Left: sketch of a the prismatic rough elements, with h the element's height, its width.Middle : top view of the periodic pattern for the rough bed.L is the side length of the pattern's square surface area.Right: sketch showing the locations of the PIV laser sheets.Positions of the green continuous lines were used in all configurations, with local coordinates y = /2, 3/4 , , 1/4(L − ) and 1/2(L − ) (from the top line to the bottom line).Positions of the black dash-lines were used in only one configuration (k/ = 1), with local coordinates y = 5/8 , 7/8 , 1/8(L − 3 ), 3/8(L − ) (from the top line to the bottom line), in order to investigate the spatial convergence of double-averaged quantities.

Figure 2 .
Figure 2. Normalized spatial convergence errors for double-averaged velocities u and w , spatiallyaveraged Reynolds shear stress u w and dispersive shear stresses ũ w for a modified-relativesubmergence β = 3.1 and an aspect ratio k/ = 1.Black diamonds ( ) are for five PIV slices and white circles (•) for three PIV slices.
sin γ and modeled as in Raupach et al.[8] under the form F x = − (1/2) C d a u2 where C d is the drag coefficient and a the frontal area per canopy volume.Vertical profiles of the total drag force F x for the three canopies have been plotted on the left part of figure4.For the elongated canopy k/ = 6, a constant non-dimensional force close to unity appears in the lower part of the canopy, corresponding to the constant velocity discussed earlier.The same region is observed with the intermediate canopy with k/ = 3.For the shallow canopy composed of cubes with k/ = 1, the non-constant distribution of drag associated with the mixing layer near the top of the rough elements reaches the bottom of the flume, and no constant velocity/drag region exists anymore.Profiles of the drag coefficient C d have been plotted on the right part of figure4for the three canopies.For k/ = 6 and 3, C d is almost constant in the lower part of the canopy, with values close to unity.In the upper parts of the canopy (i.e., (k − z)/ < 1), C d values are evolving differently for the three canopies, with values away from unity, in the range [0.5, 2].