Vegetated Channel Flows : Turbulence Anisotropy at Flow – Rigid Canopy Interface

This laboratory study aimed at investigating the mean and turbulent characteristics of a densely vegetated flow by testing four different submergence ratios. The channel bed was covered by a uniform array of aligned metallic cylinders modeling rigid submerged vegetation. Instantaneous velocities, acquired with a three-component acoustic Doppler velocimeter (ADV), were used to analyze the mean and turbulent flow structure. The heterogeneity of the flow field was described by the distributions of mean velocities, turbulent intensities, skewness, kurtosis, Reynolds stresses, and Eulerian integral scales. The exchange processes at the flow–vegetation interface were explored by applying the turbulence triangle technique, a far less common technique for vegetated flows based on the invariant maps of the anisotropic Reynolds stress tensor.


Introduction
Aquatic and riparian vegetation is a fundamental feature of riverine systems, and its presence deeply affects the hydraulic resistance, mass and momentum transport, and turbulent flow structure [1,2] with a variety of ecological and physical implications [3,4].The presence of a submerged macro-roughness represented by a vegetative obstruction deeply alters the mean and turbulent structure of the flow with implications for mass and momentum exchange across the flow-vegetation interface.The hydrodynamic structure of the flow over a submerged canopy was investigated with reference to two main vegetation models: (1) rigid vegetation, in which the vegetative obstruction is generally represented by an array of rigid cylinders or prismatic elements [5][6][7], and (2) flexible vegetation, generally simulated by blade-shaped vegetation [8,9].
The rigid cylinder model is not able to reproduce a variety of flow-influencing mechanisms exhibited by natural vegetation, generally presenting a wide range of different biomechanical properties [10] and velocity-dependent drag characteristics [11,12].Nevertheless, rigid cylinder models can be considered qualitatively effective for simulating the cross-section scale effects on the turbulence structure and describing morphodynamic processes of channels with vegetated floodplains [13].Because of the practical importance of modeling turbulence in natural watercourses and coastal wetlands, an increasing interest in describing the turbulent flow structure in vegetated flows [14][15][16], highly rough beds, and obstructed flows [17] is noticeable.
The flow over rigid submerged vegetation has been widely investigated in order to assess the effects of vegetation (density and submergence) on the flow structure and the implications for the hydraulic resistance, turbulent structures, mixing processes, and sediment transport [5,6,[18][19][20][21].A distinctive feature of such vegetated shear flows is represented by the anisotropy of Reynolds stresses, Geosciences 2018, 8, 259; doi:10.3390/geosciences8070259www.mdpi.com/journal/geoscienceswhich reflects, to a large extent, the presence of organized motions.Indeed, the anisotropic Reynolds stress component is the only part of the total stress responsible for the momentum transport [22], and, therefore, its study is crucial for the understanding and modeling of turbulence in vegetated flows.The anisotropy pattern within the flow domain can be investigated using the technique of the anisotropy invariants proposed by Lumley and Newman [23].This technique, fairly less common for the analysis of vegetated flows, allows for the characterization of the spatial distribution of the anisotropy degree and nature.Only few examples of the application of this methodology are available in the literature, and these mainly refer to numerical and experimental analyses of flows over highly rough beds [24,25].
In this study, the turbulent structure of the flow over a submerged array of rigid cylinders was experimentally investigated by coupling the traditional approach based on the spatial distributions of velocity statistics and spectral analysis [5,6,26] with a new methodology adopted for providing the overall description of the turbulence anisotropy, represented by the turbulence triangles [27].Specifically, the effects of the variability of the submergence, that is, the ratio between the uniform flow depth and the cylinder's height, on the mean velocity profiles, the distributions of higher-order velocity statistics, quadrant analysis, and the distribution of integral time and length scales were analyzed with reference to the anisotropy pattern resulting from the invariant maps.

Laboratory Flume
The experiments were carried out in an 8 m long, 0.4 m wide, and 0.4 m high Plexiglas-walled recirculating flume in the Laboratory of Hydraulics and Hydraulic Structures at the University of Naples Federico II (Figure 1a).The channel was supplied by a 4.5 m 3 tank into which water was led through the water supply system of the laboratory from a water tower.Specifically, in order to stabilize flow rates in the channel, water was not directly pumped into the flume but, from the recirculating system reservoir, was led to a water tower with a pumping station.
The discharge was measured by a magnetic flow meter (to an accuracy of ±0.1 L/s) installed on the feeding pipeline of the channel.The discharge ranged up to 30 L/s.The channel slope was set to 1%.At the inlet section, a parabolic transition from the feeding tank to the experimental channel was placed in order to reduce the disturbance at the inlet and to damp the related turbulence.which reflects, to a large extent, the presence of organized motions.Indeed, the anisotropic Reynolds stress component is the only part of the total stress responsible for the momentum transport [22], and, therefore, its study is crucial for the understanding and modeling of turbulence in vegetated flows.The anisotropy pattern within the flow domain can be investigated using the technique of the anisotropy invariants proposed by Lumley and Newman [23].This technique, fairly less common for the analysis of vegetated flows, allows for the characterization of the spatial distribution of the anisotropy degree and nature.Only few examples of the application of this methodology are available in the literature, and these mainly refer to numerical and experimental analyses of flows over highly rough beds [24,25].
In this study, the turbulent structure of the flow over a submerged array of rigid cylinders was experimentally investigated by coupling the traditional approach based on the spatial distributions of velocity statistics and spectral analysis [5,6,26] with a new methodology adopted for providing the overall description of the turbulence anisotropy, represented by the turbulence triangles [27].Specifically, the effects of the variability of the submergence, that is, the ratio between the uniform flow depth and the cylinder's height, on the mean velocity profiles, the distributions of higher-order velocity statistics, quadrant analysis, and the distribution of integral time and length scales were analyzed with reference to the anisotropy pattern resulting from the invariant maps.

Laboratory Flume
The experiments were carried out in an 8 m long, 0.4 m wide, and 0.4 m high Plexiglas-walled recirculating flume in the Laboratory of Hydraulics and Hydraulic Structures at the University of Naples Federico II (Figure 1a).The channel was supplied by a 4.5 m 3 tank into which water was led through the water supply system of the laboratory from a water tower.Specifically, in order to stabilize flow rates in the channel, water was not directly pumped into the flume but, from the recirculating system reservoir, was led to a water tower with a pumping station.
The discharge was measured by a magnetic flow meter (to an accuracy of ±0.1 L/s) installed on the feeding pipeline of the channel.The discharge ranged up to 30 L/s.The channel slope was set to 1%.At the inlet section, a parabolic transition from the feeding tank to the experimental channel was placed in order to reduce the disturbance at the inlet and to damp the related turbulence.All the tests were carried out under uniform flow conditions with different cylinder submergences.In the range of considered flow rates, a uniform flow reach was identified by analyzing the backwater profiles, experimentally determined using a gauging needle (with an accuracy of 0.1 mm).For all the test runs, uniform flow conditions were observed to be restored at All the tests were carried out under uniform flow conditions with different cylinder submergences.In the range of considered flow rates, a uniform flow reach was identified by analyzing the backwater profiles, experimentally determined using a gauging needle (with an accuracy of 0.1 mm).For all the test runs, uniform flow conditions were observed to be restored at approximately 4.5 m from the outlet section.Measurements were taken at 2.5 m from the channel inlet, within the uniform flow section of the flume and sufficiently downstream from the inlet.

Rigid Vegetation Model
This experimental study was performed by considering a uniformly vegetated channel bed in which metallic cylinders were used to simulate rigid, submerged vegetation (Figure 1b).The 45 mm long rods of 4 mm diameter, installed in holes bored into the Plexiglas bed panels, were arranged in an aligned pattern with a constant density all over the channel bed.The spacing between cylinders, in the streamwise L x and spanwise L y directions, was equal to 25 mm (Figure 1b).Consequently, the number of cylinders per unit bed area n was equal to 1600 m −2 .The vegetation density was described by three other parameters: • the frontal area per canopy volume a = d/ L x L y , equal to 6.4 m −1 ; • the frontal area per bed area λ, also known as the density roughness, computable as ak for vertically uniform vegetation, and equal to 0.288; • the solid volume fraction occupied by the canopy elements φ, evaluable as (π/4)ad for cylindrical elements and equal to 0.020.
The considered vegetation geometry and density, typical of dense canopies [28], were comparable to those commonly exhibited by aquatic submerged vegetation, as reported in Table 1.

Test Cases and Velocity Measurements
The effects of the submergence ratio on the flow structure were experimentally investigated for the hydraulic conditions reported in Table 2. Specifically, four different submergence ratios were considered: 2.4, 2.8, 3.1, and 3.4.For each condition, the canopy edge velocity (U k ), the measured friction velocity at the canopy edge (U * = −u w ), and the bulk velocity (U b = Q/Bh) were evaluated.Furthermore, the relevant Reynolds numbers (the roughness Re k = kU * /ν, flow Re h = hU b /ν, and vegetation element Re d = dU * /ν), and the flow Froude number are indicated (Table 2).Under the considered conditions of shallow submerged canopies [28,29], with h/k < 5, both turbulent stress and potential gradients drive the flow within the canopy [30].
Instantaneous velocity measurements were performed with a SonTek 16 MHz MicroADV (SonTek, San Diego, CA, USA), a three-axis acoustic Doppler velocimeter (ADV) with a down-looking probe.In order to describe the turbulence structure within and above the canopy, the maximum available sampling frequency of the instrument was adopted, that is, 50 Hz.In order to achieve well-converging higher-order moments, a sample size of more than 50,000 values was assumed [31].
ADV data were filtered by assuming a threshold correlation and signal-to-noise-ratio level of 70% and 30 dB, respectively.Finally, data were de-spiked using the phase-space threshold method [32] as modified by Wahl [33].Velocity measurements were performed in the uniform flow reach, along a vertical line in the channel midline (Figure 2a,b).The sampling volume (SV) was approximately cylindrical with a 4 mm diameter and 4.5 mm height (Figure 2a); therefore, the distances between each acquisition point were selected to be equal to 10 mm, in order to avoid overlapping the SVs of each measurement point.( / ) ( ) ADV data were filtered by assuming a threshold correlation and signal-to-noise-ratio level of 70% and 30 dB, respectively.Finally, data were de-spiked using the phase-space threshold method [32] as modified by Wahl [33].Velocity measurements were performed in the uniform flow reach, along a vertical line in the channel midline (Figure 2a,b).The sampling volume (SV) was approximately cylindrical with a 4 mm diameter and 4.5 mm height (Figure 2a); therefore, the distances between each acquisition point were selected to be equal to 10 mm, in order to avoid overlapping the SVs of each measurement point.Considering the distance between the center of the SV and the acoustic transmitter of the probe (50 mm), the flow was investigated up to 50 mm from the water surface (Figure 2a).The orientation of the probe was defined to identify the streamwise axis of the channel with the x-axis of the probe; y is the spanwise direction (positive leftwards), and z is the vertical direction (positive upwards).The good alignment and the verticality of the instrument were addressed with a laser cross-level and a three-dimensional (3D) bubble level.

Mean Velocity Profiles
In Figure 3, the mean velocity profiles of , , and (for , , and directions, respectively) are shown.The velocity and distance from the bed were made dimensionless with the streamwise mean velocity at the canopy top and the vegetation height, respectively.Considering the distance between the center of the SV and the acoustic transmitter of the probe (50 mm), the flow was investigated up to 50 mm from the water surface (Figure 2a).The orientation of the probe was defined to identify the streamwise axis of the channel with the x-axis of the probe; y is the spanwise direction (positive leftwards), and z is the vertical direction (positive upwards).The good alignment and the verticality of the instrument were addressed with a laser cross-level and a three-dimensional (3D) bubble level.

Mean Velocity Profiles
In Figure 3, the mean velocity profiles of U, V, and W (for x, y, and z directions, respectively) are shown.The velocity and distance from the bed were made dimensionless with the streamwise mean velocity at the canopy top and the vegetation height, respectively.The selected scaling velocity parameters and the resulting velocity profiles shown in Figure 3 were consistent with the results of previous studies [5,6,34,35].The discontinuity in drag at the vegetation top modified the longitudinal velocity profile, characterized by an inflection point at the top of the canopy, resembling a canonical mixing layer [28].Specifically, the drag exerted by the canopy decelerated the flow in the vegetated layer; as firstly addressed by Raupach et al. [34], this effect dominates the transfer of momentum across the flow-vegetation interface.This effect, slightly more pronounced for higher values of ℎ/ , is characteristic of dense canopies and, specifically, The selected scaling velocity parameters and the resulting velocity profiles shown in Figure 3 were consistent with the results of previous studies [5,6,34,35].The discontinuity in drag at the vegetation top modified the longitudinal velocity profile, characterized by an inflection point at the top of the canopy, resembling a canonical mixing layer [28].Specifically, the drag exerted by the canopy decelerated the flow in the vegetated layer; as firstly addressed by Raupach et al. [34], this effect dominates the transfer of momentum across the flow-vegetation interface.This effect, slightly more pronounced for higher values of h/k, is characteristic of dense canopies and, specifically, canopies with ak >≈ 0.23 [28].The large interfacial coherent structures, also called canopy-scale vortices, occupy the surface layer and partly penetrate the vegetation layer, depending on the canopy density.Nepf et al. [35] defined a parameter describing the scaling for the penetration depth (δ e ): where C D is the drag coefficient of vegetation elements.In this study, for C D = 1, δ e varied between 27 and 45 mm and δ e /k varied between 0.6 and 1.
The turbulent flow in the investigated region of the channel was fully developed, and, considering that the leading edge of the canopy was located approximately 2.5 m upstream from the measurement section, the canopy-scale vortices and the penetration depth reached a stable condition.

Turbulent Intensities
In Figure 4, the turbulent intensities in the three directions, defined as the root-mean-square values for the three velocity components, made dimensionless with the measured friction velocity at the canopy top, are shown.
The profiles showed a similar trend, collapsing together along the same curve.The maximum longitudinal turbulent intensities occurred at the top and directly above the canopy.The turbulent intensity decreased within the vegetation layer and towards the free surface.Indeed, the vegetation damped the turbulence, and this effect was observed to be slightly more pronounced for higher submergence ratios.At the canopy top, the normalized vertical fluctuating velocity component was ≈1 for each submergence.For dense canopies, the characteristic value is about 1.1, which is typical for mixing layers, in agreement with Raupach et al. [34] and Poggi et al. [5].Larger values, up to 1.3, are typical for sparse canopies.At the canopy top, the normalized streamwise and vertical fluctuating velocity components attained values in agreement with Raupach et al. [34], Poggi et al. [5], and Nezu et al. [6] for dense canopies.

Skewness and Kurtosis
The skewness of the velocity fluctuations is shown in Figure 5.A non-zero skewness indicates an asymmetric probability density function of the considered variable, which means that larger excursions in one direction were more probable than in the other direction, depending on the sign of the statistics.The profiles of Figure 5

Skewness and Kurtosis
The skewness of the velocity fluctuations is shown in Figure 5.A non-zero skewness indicates an asymmetric probability density function of the considered variable, which means that larger excursions in one direction were more probable than in the other direction, depending on the sign of the statistics.The profiles of Figure 5 describe a clear trend for the skewness in u to take positive values inside the canopy and negative values outside and, complementarily, for the skewness in w to take negative values inside the canopy and positive values outside.This trend indicates the dominant role of sweep events inside the cylinder array and the ejection events above the canopy, in agreement with the results of other experimental analyses [5,6,14].The intensity of these effects, as confirmed by quadrant analysis, increased with the submergence ratio because of the increasing momentum transfer between the vegetated and non-vegetated zones.A normal distribution could be assumed for the fluctuations in the spanwise direction.

Skewness and Kurtosis
The skewness of the velocity fluctuations is shown in Figure 5.A non-zero skewness indicates an asymmetric probability density function of the considered variable, which means that larger excursions in one direction were more probable than in the other direction, depending on the sign of the statistics.The profiles of Figure 5 describe a clear trend for the skewness in to take positive values inside the canopy and negative values outside and, complementarily, for the skewness in to take negative values inside the canopy and positive values outside.This trend indicates the dominant role of sweep events inside the cylinder array and the ejection events above the canopy, in agreement with the results of other experimental analyses [5,6,14].The intensity of these effects, as confirmed by quadrant analysis, increased with the submergence ratio because of the increasing momentum transfer between the vegetated and non-vegetated zones.A normal distribution could be assumed for the fluctuations in the spanwise direction.The kurtosis profiles evaluated for the three velocity fluctuation components, shown in Figure 6, consistent with the results of Poggi et al. [5] for dense canopies, tended to assume higher values in the regions immediately around the canopy top.This effect, as confirmed by the quadrant z/k The kurtosis profiles evaluated for the three velocity fluctuation components, shown in Figure 6, consistent with the results of Poggi et al. [5] for dense canopies, tended to assume higher values in the regions immediately around the canopy top.This effect, as confirmed by the quadrant analysis illustrated in the next section, confirmed the dominant role of sweep and ejection events for the Reynolds shear-stress production and the momentum transport at the canopy top.Within the canopy, was observed to attain values of >1, consistent with the lateral fluctuations related to the von Karman vortex shedding process induced by the presence of cylinders.

Reynolds Shear Stress
In Figure 7 Within the canopy, Ku v was observed to attain values of >1, consistent with the lateral fluctuations related to the von Karman vortex shedding process induced by the presence of cylinders.

Reynolds Shear Stress
In Figure 7, the Reynolds shear stresses −u i u j , normalized by the squared friction velocity at the canopy top, are plotted along the inspected vertical for the four different submergences.The Reynolds stress −u w profile attained a peak at the canopy top (z/k = 1) and exhibited a sharp decrease within the canopy.This effect, as observed by other authors [6,30], was due to the presence of the canopy elements, which inhibited the momentum transfer between the surface layer and the underlying vegetated layer.Within the canopy, was observed to attain values of >1, consistent with the lateral fluctuations related to the von Karman vortex shedding process induced by the presence of cylinders.

Reynolds Shear Stress
In Figure 7, the Reynolds shear stresses − , normalized by the squared friction velocity at the canopy top, are plotted along the inspected vertical for the four different submergences.The Reynolds stress − ′ ′ profile attained a peak at the canopy top ( / = 1) and exhibited a sharp decrease within the canopy.This effect, as observed by other authors [6,30], was due to the presence of the canopy elements, which inhibited the momentum transfer between the surface layer and the underlying vegetated layer.The penetration depth (Table 3), estimated as the elevation at which the Reynolds stress − ′ ′ decays to 10% of its maximum value [30], was evaluated from Figure 7 for the different submergence z/k The penetration depth (Table 3), estimated as the elevation at which the Reynolds stress −u w decays to 10% of its maximum value [30], was evaluated from Figure 7 for the different submergence ratios, showing a good agreement with the estimation of the scaling of Equation 1. Depending only on the vegetation density, approximately constant values were observed for all the test runs.As the submergence increased, the shape of the normalized xy shear-stress profile changed, exhibiting higher values.These effects were ascribable to the growing effects of the secondary currents inside the experimental channel [6], considering that the aspect ratio of the channel (B/(h − k)) ranged from 3.8 to 6.2.

Quadrant Analysis
The trend emerging from the skewness and kurtosis distributions was confirmed by the results of the quadrant analysis.This simple conditional-sampling technique can provide useful information for the interpretation and the detection of organized coherent structures in the flow, allowing the nature of contributions to the Reynolds stress u w from the different events to be understood [36].In quadrant analysis [36][37][38], the velocity fluctuation components are plotted on a u − w plane.The plane is divided into four quadrants; each one is characteristic of a different event: • quadrant 1: outward interaction (u > 0, w > 0); • quadrant 2: ejection (u < 0, w > 0); • quadrant 3: inward interaction (u < 0, w < 0); • quadrant 4: sweep (u > 0, w < 0).
In order to highlight the contribution of extreme events to the overall Reynolds stress, a hyperbolic threshold was introduced, given by the expression |u w | = Hu rms w rms [5,6,36], which defines, on the u − w plane, a hyperbolic hole; the parameter H, the size of the hole, was assumed to be equal to 2 [36].At any point in the flow domain, the total contribution of each quadrant to the Reynolds stress can be calculated as follows [18]: where the subscript i refers to the ith quadrant and I i is a dummy variable equal to 1 if u (t)w (t) belongs to the ith quadrant and equal to 0 otherwise.Clearly, The contribution of the extreme events of each quadrant, instead, can be calculated as follows [36]: where H is the threshold level defining the size of the hyperbolic hole region and the subscript i refers to the ith quadrant.I i,H is a dummy variable equal to 1 if u (t)w (t) belongs to the ith quadrant and |u (t)w (t)| > Hu rms w rms and equal to 0 otherwise.In Figure 8, the contribution of each quadrant to the Reynolds stress is plotted as a function of the normalized distance from the bed.The contributions of the different events were firstly evaluated considering H equal to 0 and secondly by considering H equal to 2 (Figure 8), to make more evident the role played by sweep and ejection events in the Reynolds stress production.The quadrant analysis confirmed the central role of sweep within the canopy and, complementarily, the dominant role of ejection events in the region above the canopy, only slightly affected by the submergence.Inward and outward interaction event contributions were almost equal and negligible.At / ≈ 1.5, the contributions of sweep and ejection were approximately equal.The quadrant analysis confirmed the central role of sweep within the canopy and, complementarily, the dominant role of ejection events in the region above the canopy, only slightly affected by the submergence.Inward and outward interaction event contributions were almost equal and negligible.At z/k ≈ 1.5, the contributions of sweep and ejection were approximately equal.

Turbulence Anisotropy
In order to provide an overall description of turbulence in the flow over the rigid vegetated bed, mean time invariant analyses were performed using the turbulence triangle [23,27], an anisotropy-invariant map.The anisotropy of turbulence is described by the Reynolds stress anisotropy tensor: The non-dimensional form of the anisotropy tensor, recognized in the trace of the Reynolds stress tensor as twice the turbulence kinetic energy (K = 1 2 u i u i ), is given by The invariants of b ij can quantify the different states of turbulence and the relative strength of the fluctuating velocity components, that is, the componentality of the turbulence field [39].
Anisotropy-invariant maps are two-dimensional (2D) domains based on invariant properties of b ij .Specifically, in the turbulence triangles, the axis variables are linear combinations of the Reynolds stress anisotropy eigenvalues λ i ; the map defines the domain of all possible turbulent flows.Each boundary of the map is characteristic of a particular turbulence state (vertices are related to one-dimensional (1D), 2D, and 3D turbulence; edges are related to specific turbulent processes: axisymmetric expansion, axisymmetric contraction, and two-component turbulence).More details on linear and non-linear invariable maps and the shape and meaning of boundaries can be found in [40].
In Figure 9, the structure of the flow over the rigid vegetated bed is described through a turbulence triangle.Each point has coordinates η = √ I I/3 and ξ = 3 √ I I I/2 [40], where I I and I I I are the quadratic and cubic invariants of the Reynolds stress anisotropic tensor, respectively.
The patterns defined by the points of the invariant map, in the range of the considered Reynolds numbers and submergences, were clear and approximately equal for all the considered conditions.
Moving from the bed toward the water surface, the points of the invariants of the Reynolds stress tensor described a peculiar path.Near the channel bed, at z/k approximately equal to 0.1, the turbulence approached a 2D state.Moving upwards, the turbulence reached, at the top of the canopy (z/k = 1), a quasi-1D state: the fluctuation in the streamwise direction was dominant because of the presence of the large-scale coherent structures, and the shape of the turbulence was cigar-like.Specifically, for all the investigated conditions, the 1D characteristic and the degree of anisotropy were maximums at z/k = 1.Above the canopy (z/k > 1), the path described by the points of the invariants moved through the axisymmetric expansion tending to an isotropic state.Because of the dimension of the ADV SV, for each submergence, the last measurement point was 5 cm from the surface (for the highest flow condition, the last point was at z/k equal to 2.1; the water surface was at z/k equal to 3.4); therefore, the evolution of the path to the water surface is not shown.Nevertheless, because of the presence of the water surface, which inhibits the vertical fluctuations, the turbulence structure was expected to become 2D [25].
The 2D state of the turbulence inside the canopy, and, specifically, in the lower part of the vegetated layer, indicated the growing importance of the lateral fluctuations due to the von Karman vortex shedding process triggered by the cylinder array, as shown also by the plots of Figure 6.The onset of this process was ensured by the local high stem Reynolds number Ud/υ ≈ 850 > 150-200 [41].
In Figure 10, the dimensionless power density spectra of the lateral velocity fluctuations, evaluated at z/k = 0.1 for all the investigated conditions, are shown together.The spectral density was normalized by dividing by the product of the root-mean-square level of the local longitudinal velocity and the stem diameter.The frequency is given in terms of Strouhal number (St = f d/U).A concentration of energy in the range of the Strouhal number between 0.2 and 0.1 was observed, in agreement with the frequency of the vortex shedding process observed by Poggi et al. [5] and Zong et al. [42] for cylinder arrays.The von Karman vortex shedding peak in the power density spectra was only slightly visible because of the deep penetration of the canopy-scale structures within the array (Table 3).The patterns defined by the points of the invariant map, in the range of the considered Reynolds numbers and submergences, were clear and approximately equal for all the considered conditions.
Moving from the bed toward the water surface, the points of the invariants of the Reynolds stress tensor described a peculiar path.Near the channel bed, at / approximately equal to 0.1, the turbulence approached a 2D state.Moving upwards, the turbulence reached, at the top of the canopy ( / = 1), a quasi-1D state: the fluctuation in the streamwise direction was dominant because of the presence of the large-scale coherent structures, and the shape of the turbulence was cigar-like.Specifically, for all the investigated conditions, the 1D characteristic and the degree of anisotropy were maximums at z/k = 1.Above the canopy ( / > 1), the path described by the points of the invariants moved through the axisymmetric expansion tending to an isotropic state.Because of the dimension of the ADV SV, for each submergence, the last measurement point was 5 cm from the surface (for the highest flow condition, the last point was at / equal to 2.1; the water surface was at / equal to 3.4); therefore, the evolution of the path to the water surface is not shown.Nevertheless, because of the presence of the water surface, which inhibits the vertical fluctuations, the turbulence structure was expected to become 2D [25].
The 2D state of the turbulence inside the canopy, and, specifically, in the lower part of the vegetated layer, indicated the growing importance of the lateral fluctuations due to the von Karman vortex shedding process triggered by the cylinder array, as shown also by the plots of Figure 6.The onset of this process was ensured by the local high stem Reynolds number / ≈ 850 > 150-200 [41].
In Figure 10, the dimensionless power density spectra of the lateral velocity fluctuations, evaluated at / = 0.1 for all the investigated conditions, are shown together.The spectral density was normalized by dividing by the product of the root-mean-square level of the local longitudinal  3).In order to highlight the effects of the rigid cylinders on the turbulent structure, a comparison between the anisotropy-invariant maps of the vegetated bed (i.e., Figure 9d) and the correspondent maps from LES applied to an open channel flow over a smooth bed at = 13,680 [25] is shown in Figure 11.In order to highlight the effects of the rigid cylinders on the turbulent structure, a comparison between the anisotropy-invariant maps of the vegetated bed (i.e., Figure 9d) and the correspondent maps from LES applied to an open channel flow over a smooth bed at Re = 13,680 [25] is shown in Figure 11.In order to highlight the effects of the rigid cylinders on the turbulent structure, a comparison between the anisotropy-invariant maps of the vegetated bed (i.e., Figure 9d) and the correspondent maps from LES applied to an open channel flow over a smooth bed at = 13,680 [25] is shown in Figure 11.Although the two plots refer to flows with different Reynolds numbers (13,680 and 75,000 for the smooth and uniformly vegetated beds, respectively), the comparison allows the main differences in terms of the turbulence anisotropy pattern to be seen.While in the viscous sublayer, directly above the flume bed, a two-component isotropic turbulence was expected for both the conditions; at / = 0.1 ( /ℎ ≈ 0.04) the behavior was different and the turbulence in the vegetated case presented a 2D structure.On the contrary, in the smooth bed, as expected for a canonical boundary layer, the Although the two plots refer to flows with different Reynolds numbers (13,680 and 75,000 for the smooth and uniformly vegetated beds, respectively), the comparison allows the main differences in terms of the turbulence anisotropy pattern to be seen.While in the viscous sublayer, directly above the flume bed, a two-component isotropic turbulence was expected for both the conditions; at z/k = 0.1 (z/h ≈ 0.04) the behavior was different and the turbulence in the vegetated case presented a 2D structure.On the contrary, in the smooth bed, as expected for a canonical boundary layer, the turbulence was highly anisotropic and 1D.Moving upwards, the evolution of the pattern was still different, and, while for the vegetated bed the anisotropy progressively increased and a 1D characteristic was observed, for the smooth bed, the turbulence gradually returned to isotropy via an axisymmetric expansion process.For the vegetated bed, turbulence at z/k ≈ 2 (z/h ≈ 0.6) still presented a high degree of anisotropy and a cigar-like shape, because this area of the flow was dominated by the canopy-scale coherent structures.For the smooth bed, the turbulence became quasi-isotropic at z/h ≈ 0.7 and then presented, approaching the free surface, a tendency toward a 2D structure.

Integral Scales
In order to investigate the time and spatial extent of the coherent structures dominating the flow, the calculation of Eulerian integral scales was performed.Specifically, an integral time scale is a measure of the memory effect in the flow field of the persistence of a large-scale eddy at a fixed Eulerian point.In Figure 12a, integral time scales evaluated by means of Equation ( 6) for the three directions along the inspected vertical and for the different considered conditions are shown.
The integration domain for the determination of the integral time scale ranged from zero up to the value at which the autocorrelation function fell to 1/e [43,44] so as to define an integration domain independent of the shape of the autocorrelation function [45].Large eddies persist for a long time and are advected slowly; thus the integral scale was large.
The Eulerian integral length scales could be estimated from the single-point integral time scales by applying Taylor's frozen-turbulence hypothesis and assuming the mean local longitudinal velocity U as the convection velocity U c of the mean eddies using the equation L E,i = UT E,i .
The integral length scales evaluated for the four different conditions are shown in Figure 12b.The results confirmed the interpretation of the flow structure drawn from the anisotropy analysis.
value at which the autocorrelation function fell to 1/e [43,44] so as to define an integration domain independent of the shape of the autocorrelation function [45].Large eddies persist for a long time and are advected slowly; thus the integral scale was large.The Eulerian integral length scales could be estimated from the single-point integral time scales by applying Taylor's frozen-turbulence hypothesis and assuming the mean local longitudinal velocity as the convection velocity of the mean eddies using the equation , = , .The integral length scales evaluated for the four different conditions are shown in Figure 12b.The results confirmed the interpretation of the flow structure drawn from the anisotropy analysis.The turbulence exhibited a quasi-2D structure near the channel bed, within the vegetation layer, as a result of the dominant effects of the stem-scale turbulence and von Karman vortices.Moving toward the vegetation top, while , tended to decrease, , and , showed a progressively increasing trend, confirming the axisymmetric expansion process shown by the invariant maps.Specifically, at the canopy top, , and , were approximately equal to and 0.4 , showing that, in the upper vegetated layer, the dominant eddy size was in the order of , consistent with the results of Raupach et al. [34] and Nezu et al. [6].Over the canopy, the three integral length scales became larger; in this region the flow was dominated by the canopy-scale turbulence.
In Figure 13, the good fit between the integral length scales relative to runs 3 and 4 and the results of Nezu et al. for dense canopies [6] are shown.The turbulence exhibited a quasi-2D structure near the channel bed, within the vegetation layer, as a result of the dominant effects of the stem-scale turbulence and von Karman vortices.Moving toward the vegetation top, while L E,y tended to decrease, L E,x and L E,z showed a progressively increasing trend, confirming the axisymmetric expansion process shown by the invariant maps.Specifically, at the canopy top, L E,x and L E,z were approximately equal to k and 0.4k, showing that, in the upper vegetated layer, the dominant eddy size was in the order of k, consistent with the results of Raupach et al. [34] and Nezu et al. [6].Over the canopy, the three integral length scales became larger; in this region the flow was dominated by the canopy-scale turbulence.
In Figure 13, the good fit between the integral length scales relative to runs 3 and 4 and the results of Nezu et al. for dense canopies [6] are shown.
that, in the upper vegetated layer, the dominant eddy size was in the order of , consistent with the results of Raupach et al. [34] and Nezu et al. [6].Over the canopy, the three integral length scales became larger; in this region the flow was dominated by the canopy-scale turbulence.
In Figure 13, the good fit between the integral length scales relative to runs 3 and 4 and the results of Nezu et al. for dense canopies [6] are shown.In comparison with terrestrial canopy, aquatic canopy is depressed by the free surface, while, on the contrary, terrestrial canopy is unconfined.This result confirms how, in aquatic canopy, eddies may be influenced by the submergence.

Conclusions
The results of this experimental analysis provided an exhaustive picture of the effects of increasing submergence on the mean and turbulent structure of the flow over an array of submerged rigid cylinders modeling dense vegetation.Specifically, rigid cylinders were arranged in an aligned pattern, modeling submerged rigid vegetation.Different submergence ratios were investigated.
A variation of approximately 50% of h/k was observed not to significantly affect the flow structure, indicating that, in the range of the tested submergences, the stem density played a major role in the hydrodynamic structure of the flow, as reported by Nezu et al. and Poggi et al. [5,6].
The mean velocity profiles presented a characteristic inflected shape because of the different flow velocities within and above the vegetation; this was slightly more pronounced for the higher values of submergence.This effect, related to the high density of the cylinder array, was the main feature of the considered obstructed flow, and the consequent Kelvin-Helmholtz-type instability dominated the momentum exchange across the vegetated interface.
The higher-order moments and the quadrant analysis allowed the Reynolds stress production process and distribution throughout the water depth to be characterized, confirming the dominant role of sweep and ejection bursting events.The maximum Reynolds stress was observed at the top of the canopy.A sudden decay of τ zx within the vegetation, due to the significant contribution of the drag force to the momentum balance, was observed.The integral length scales showed a trend characteristic of dense aquatic canopy, allowing eddies of the flow field to be characterized.
The invariant maps, generally used for verification of numerical model results and here applied for the interpretation of the componentality of the turbulence in the vegetated flow, gave very interesting results, showing their usefulness in interpreting the turbulent structure of vegetated flows.In particular, the flow structure picture emerging from the anisotropy analysis was consistent with the traditional statistical analysis of turbulence, providing a complete description of the flow field.Specifically, the method can be easily transferred to more complex vegetation models and configurations and, together with traditional analysis, can contribute to the understanding of turbulence in vegetated contexts.

Figure 1 .
Figure 1.(a) Experimental facility.The laboratory flume was 8 m long with a 0.4 m wide and 0.4 m deep rectangular cross-section.The bed of the flume was uniformly covered by an array of rigid cylinders.(b) Geometry of the array of aligned cylinders and specification of diameter (d), spacing ( and ), and cylinder height (k).

Figure 1 .
Figure 1.(a) Experimental facility.The laboratory flume was 8 m long with a 0.4 m wide and 0.4 m deep rectangular cross-section.The bed of the flume was uniformly covered by an array of rigid cylinders.(b) Geometry of the array of aligned cylinders and specification of diameter (d), spacing (L x and L y ), and cylinder height (k).

Figure 2 .
Figure 2. (a) Position of measurement stations and geometry of sampling volume.(b) Location of inspected vertical within the cylinder array.

Figure 2 .
Figure 2. (a) Position of measurement stations and geometry of sampling volume.(b) Location of inspected vertical within the cylinder array.

16 Figure 3 .
Figure 3. Dimensionless mean velocity profiles for U, V, and W. Dash-dotted line indicates the canopy top.

Figure 3 .
Figure 3. Dimensionless mean velocity profiles for U, V, and W. Dash-dotted line indicates the canopy top.

Figure 4 .
Figure 4. Turbulent intensities in the three directions x, y, and z.

Figure 4 .
Figure 4. Turbulent intensities in the three directions x, y, and z.

Figure 4 .
Figure 4. Turbulent intensities in the three directions x, y, and z.

Figure 5 .
Figure 5. Skewness of the fluctuation velocity components in the three directions x, y, and z.

Figure 5 .
Figure 5. Skewness of the fluctuation velocity components in the three directions x, y, and z.
Geosciences 2018, 8, x FOR PEER REVIEW 7 of 16analysis illustrated in the next section, confirmed the dominant role of sweep and ejection events for the Reynolds shear-stress production and the momentum transport at the canopy top.

Figure 6 .
Figure 6.Kurtosis of the fluctuation velocity components in the three directions x, y, and z.

Figure 6 .
Figure 6.Kurtosis of the fluctuation velocity components in the three directions x, y, and z.

Figure 6 .
Figure 6.Kurtosis of the fluctuation velocity components in the three directions x, y, and z.

Figure 7 .
Figure 7. Reynolds stress profiles.The quantities are normalized with square of the friction velocity at the canopy top * .

Figure 7 .
Figure 7. Reynolds stress profiles.The quantities are normalized with square of the friction velocity at the canopy top U * .

Figure 8 .
Figure 8. Contribution of each quadrant to the overall Reynolds stress: (a) evaluated with Equation (2) considering all the events, and (b) evaluated with Equation (3) considering only the extreme events, assuming the parameter equal to 2. Inward and outward interactions are represented by solid squares (■) and triangles (▲), respectively, and sweep and ejection interactions by diamonds (♦) and circles (•).The different colors refer to the four different considered conditions: blue: 30 L/s; green: 25 L/s; yellow: 20 L/s; red: 15 L/s.

Figure 8 .
Figure 8. Contribution of each quadrant to the overall Reynolds stress: (a) evaluated with Equation (2) considering all the events, and (b) evaluated with Equation (3) considering only the extreme events, assuming the H parameter equal to 2. Inward S 1 and outward S 3 interactions are represented by solid squares ( ) and triangles ( ), respectively, and sweep S 2 and ejection S 4 interactions by diamonds ( ) and circles ( ).The different colors refer to the four different considered conditions: blue: 30 L/s; green: 25 L/s; yellow: 20 L/s; red: 15 L/s.

Figure 10 .
Figure 10.Power density spectra at = 5 mm for the different runs.The ordinates were made dimensionless by dividing by the local root mean square multiplied by the diameter of the rods.

Figure 10 .
Figure 10.Power density spectra at z = 5 mm for the different runs.The ordinates were made dimensionless by dividing by the local root mean square multiplied by the diameter of the rods.

Figure 10 .
Figure 10.Power density spectra at = 5 mm for the different runs.The ordinates were made dimensionless by dividing by the local root mean square multiplied by the diameter of the rods.

Figure 12 .
Figure 12.Profiles of integral scales of the flow: (a) Eulerian integral time scales; (b) Eulerian integral length scales.

Figure 12 .
Figure 12.Profiles of integral scales of the flow: (a) Eulerian integral time scales; (b) Eulerian integral length scales.

Figure 13 .
Figure 13.Comparison with analogous literature calculation of integral length scales [6,46]: (a) integral length scale profile along the vertical direction for the x component; (b) integral length scale profile along the vertical direction for the z component.In the plot, only the profiles evaluated for the two highest submergence conditions are reported.

Figure 13 .
Figure 13.Comparison with analogous literature calculation of integral length scales [6,46]: (a) integral length scale profile along the vertical direction for the x component; (b) integral length scale profile along the vertical direction for the z component.In the plot, only the profiles evaluated for the two highest submergence conditions are reported.

Table 2 .
Details of experimental conditions.

Table 2 .
Details of experimental conditions.