The experiments were conducted in a recirculating, free-surface flow flume at the Hydraulics Laboratory of the School of Soil and Water Conservation at Beijing Forestry University. The flume was 12.4 m long (L), .3 m wide (W), and .25 m deep (H ₀). The electromagnetic flowmeter system controlled the discharge. At the same time, a honeycomb-shaped straightener at the inlet and an adjustable weir gate at the end of the flume were used to set a stable and uniform flow.

The vegetation was simulated with plastic cylinders of .006 m diameter (D), and the cylinders representing short vegetation and tall vegetation were .07 m (h ₁) and .18 m (h ₂), respectively. The cylinders were embedded in five 1 m long, .3 m wide, and .5 cm thick plexiglass plates. The vegetation section was 5-m-long, and the leading edge was 2.5 m from the flume inlet (see Figure 2A). A Cartesian coordinate is established with the frontal edge of the vegetation section as the origin. The x, y, and z coordinates denote streamwise, lateral, and vertical directions, respectively. The preliminary experiment of double-layered vegetation flow was conducted, the velocities at x = 0, .5, 1.3, 2.5, 3.6, and 4.5 m were measured, and the results showed that the longitudinal velocity at x = 3.6 m and x = 4.5 m had high consistency which means the flow behind x = 3.6 m is fully developed. Therefore, the velocity was measured at x = 4.0 m, representing the velocity of the fully developed region. Several cylinders were removed in order to measure the velocity using ADV. A side-view image of the double-layer model plants and velocity measurements for run XHS1 are shown in Figure 2B.

The specific parameters of each run are given in Table 1. The experiments include two vegetation arrangements marked as X and Y, with different density ratios α = m ₁/m ₂. The water depth was set as a higher submergence state H = 14.00 cm (HS), and a lower submergence state H = 12.80 cm (LS), under which condition the short vegetation layer was submerged and the tall vegetation layer was emergent. The numbers 1-3 denoted three different slopes J, .001, .002, and .003. According to the Reynolds number Re=(U _m R)/ν and the Froude number F _r=U _m/ g H , the experimental flow is turbulent and subcritical flow, where U _m = Q/WH is the average velocity over the section, R=(H-λ ₁ h ₁-λ ₂ H)/(m ₁ h ₁ D + m ₂ HD) is the hydraulic radius of the half-submerged double-layered vegetation flow without consideration of the impact of the bed and walls (Tang et al., 2021), H is the water depth, Q is the flow discharge, g = 9.81 ms⁻² is the acceleration of the gravity, ν = 1.01 × 10⁻⁶ m²s⁻¹ is the kinematic viscosity of the water (based on the kinematic viscosity of the water under 101.325 kPa and 20°C).

Nortek Acoustic Doppler Velocimetry (ADV) was used to measure the instantaneous velocities. The ADV has one transmitting probe and four receiving probes. The transmitting probe in the center emits acoustic signals, and then receiving probes receive the acoustic frequency changes caused by the reflection or scattering of moving particles in the water to collect the data of the sampling body 5 cm in front of the probe. After the signal is processed by the system, the instantaneous velocity, signal intensity, signal-to-noise ratio (SNR), and correlation (COR) of each measuring point can be converted into a text file for output. The ADV cannot accurately measure the range of 2–3 cm below the water surface, for the upper signal-receiving probe will be emergent. The data acquisition frequency was 200 Hz, the sample volume was 5.5 mm³, and the sampling time was 40 s. Before the experiment, the ADV was fixed on a support. The tracer particles used in the experiment are hollow plexiglass beads of 20 µm diameter, and the density is .9 gcm⁻³. Figure 3 shows that there were three measuring locations in runs X representing behind the tall cylinder (BT), free stream behind the tall and short cylinders’ gap (BTS) and the short cylinder (BS), and one more for runs Y: Free stream behind two short cylinders’ gap (BSS). The measurement was done from the bed (z = .005 m) towards the free water surface with an increment of .005 m at each measuring location.

Python processed the experimental data. The longitudinal, lateral, and vertical instantaneous velocities with a signal-to-noise ratio (SNR) greater than 15 and correlation (COR) greater than 70% are regarded as valid data, and the time-average velocity of x, y, and z direction u ¯ , v ¯ , and w ¯ are obtained by time averaging the instantaneous velocities. The time-average longitudinal velocity of each measuring location was averaged horizontally to obtain the mean longitudinal velocity of the measuring cross-section.

Figure 4 shows the longitudinal velocities at different locations of runs X and Y, and the black dotted line indicates the height of short vegetation h ₁ = .07 m. The longitudinal velocity profile can be roughly divided into three zones, as mentioned above. Zone 1A is from the bed to z = .055 m (0 < z ≤ .79h ₁), which is called the constant velocity zone, the velocity almost stays constant or has an insignificant increasing tendency. Zone B is from z = .055 m to z = .07–.08 m (.79h ₁< z ≤ 1.14h ₁). The velocity reaches the maximum at the top of the short vegetation or a slightly higher place and basically remains constant at zone 2A (1.14h ₁< z ≤ 1.71h ₁). The zoning of the longitudinal velocity profile is almost the same as the previous experimental results (Liu et al., 2010; Huai et al., 2014; Tang et al., 2021).

Figure 4 shows that, except for run XHS2, the longitudinal velocity at location BT in runs X is the largest, and the one at location BTS is the smallest. The velocity at location BT in run XHS2 ranges between the value of location BS and location BTS, which is different from the regularities of the other runs, probably due to the error caused by an experimental operation, such as the slight deviation of the measurement location. Furthermore, in runs XHS1, XLS2, and XLS3, the velocities at locations BT and BTS have a percentage variance of 5%–18%. This is most likely due to the fact that the increasing of velocity component caused by the wake vortices after stem, based on the velocity model of linear superposition (Nikora et al., 2013). In runs Y, the velocity profiles at different locations are similar, and the velocity curve is not as smooth as that of runs X, especially the measured data of run YLS1, possibly due to the more complex arrangement of the vegetation cylinders of runs Y, leading to the frequent fluctuation of the velocity.

Figure 5 displays the dimensionless mean longitudinal velocity of the measuring cross-section in each run. Every dimensionless velocity shows consistent distribution characteristic of “constant-increasing-constant” in vertical. For runs X, the velocity in zone 1A (u _1A) have a magnitude of .69–.83U _m and the velocity in zone 2A u _2A = .97–1.17U _m. For runs Y, the velocity in zone 1A u _1A = .59–.76U _m, and the velocity in zone 2A u _2A = 1.04–1.35U _m. The density ratio α of the double-layer vegetation under the same water depth and slope condition dictates the magnitude of the constant velocity within it. Here, the discharge Q is also an important factor affecting the velocity within the upper layer.

We define that the ratio of the mean velocity in zone 2A to the mean velocity in zone 1A β = u 2 A ¯ / u 1 A ¯ , where u 2 A ¯ and u 1 A ¯ is the mean experimental velocity of zone 2A and zone 1A, respectively. The α of runs Y is greater than runs X, leading to a more significant difference between the vegetation resistance and velocity distribution in zones 1A and 2A, and β of runs Y is 1.59–1.80. Table 1 shows that the velocity increases with the channel slope under the condition of the same vegetation arrangement and water depth. It can also be concluded from Figure 5 that the slope does not have much effect on the velocity ratio between zones 2A and 1A, the β of runs XHS and XLS are similar as the slope changes.

The inflection point of the velocity profile in vegetated flows can be characterized by the lateral vorticity, and it can be calculated as the following equation (Singh et al., 2019): ω = ∂ w ∂ x − ∂ u ∂ z , (10) where the fluctuation of the vertical velocity ∂ w / ∂ x is much weaker than the streamwise fluctuation ∂ u / ∂ z , so the term ∂ w / ∂ x is neglected in the present calculation. The lateral vorticity profile is shown in Figure 6. In Figures6A–G, the lateral vorticity basically stays constant, being about zero in the lower zone of the short vegetation layer and upper zone of the tall vegetation layer, and there is a spike at the junction of short and tall vegetation layer. It shows that the momentum exchange in zone B is distinctly larger than in zones 1A and 2A. The results are similar to the lateral vorticity profiles of single-layered rigid vegetation flow found in Liu et al. (2008), and double-layered vegetation flow of Liu et al. (2010) and Singh et al. (2019). The vorticity distribution perfectly corresponds to the velocity distribution in Figure 4. The lateral vorticity profiles of run YLS1 show a strong fluctuation in zones 1A and 2A, and there is an exceptional value in run XHS2 at z = .115 m, possibly due to the experimental error at some measurement points.

The application of the velocity model should base on the determination of the two parameters: The length scale of the mixing layer zone B ( δ w ) and the vegetation drag coefficient (C _D).

For simplicity, the deflection points of the mixing layer are determined using the method in Truong et al. (2019), which is used to define the edge of the penetration length in the vegetated zone and non-vegetated zone. The deflection points z ₁ and z ₂ are defined as the streamwise velocity u z 1 = 1 + 5 % ∙ u 1 A ¯ and u z 2 = 1 − 5 % ∙ u 2 A ¯ , respectively.

The drag coefficient is calculated using the measurement data in zone 1A: C D = g i 1 − λ 1 1 2 D m 1 u 1 A ¯ 2 . (11)

In runs X, the velocity model almost agrees well with the experimental data (Figures 7A–F), but underestimates or overestimates the velocity of zone 2A in runs Y (Figures 7G, H).

Moreover, in order to intuitively present the deviation of the predicted velocity from the measured velocity, the relative error (RE) is defined as follow: R E = U mod − u m e a u m e a . (12) where U _mod is the predicted velocity of the model, and u _mea is the measured velocity. Figure 8 shows the RE at each measuring point of runs X along the water depth. The RE are relatively small in general, indicating that the model is a good reflection of the flow velocity profile. However, the RE in the near-bed 1 cm zone, zone B and zone 2A are larger than those in zone 1A, because the constant velocity model in zone 1A cannot catch the influence of boundary conditions precisely, and the entire velocity model is dependence on the C _D and δ w determined by the measurement data in zone 1A, leading to a higher RE in zone B and 2A.

It is mentioned in Section 4.2 that the slope does not contribute to β, and the velocity model of zone 2A applied here denotes that the uniform velocity term difference between zone 2A and zone 1A is only related to the density of two vegetation layers, i.e., U U D 2 / U U D 1 = α 1 − λ 1 / 1 − λ 2 , the calculated values of runs X and runs Y are 1.406 and 1.715 respectively, while the β of runs X is 1.35–1.44, and 1.59–1.80 of runs Y. Therefore, it shows a small deviation in the modeling of the velocity in zone 2A. The model application results also verify that the longitudinal velocity in double-layered rigid vegetation flow follows a hyperbolic tangent distribution in the mixing layer as the distribution in single-layered submerged flexible vegetation flow, which may extend to the prediction of the mixing layer velocity in multi-layered vegetation flow; for instance, the velocity in a mixing layer can be expressed as the linear superposition of a uniform velocity in the lower layer and a hyperbolic tangent profile.

The C _D of this paper was estimated from the flow condition in zone 1A, and the dependence of the C _D on the stem Reynolds number of zone 1A Re _UD1 =(U _UD1 D)/ν for runs X is illustrated in Figure 9. It shows that C _D decreases as the Re _UD1 increases and has an increasing tendency to a growing slope J. This trend is probably due to the vortices penetrating deeper into the vegetation, and the gravitational component of the water along the slope direction is more significant with an increased slope, creating a more considerable vegetation drag. The results are consistent with the experimental data of single-layered submerged vegetation flow published in Nikora et al. (2013).

This study is concerned with the influence of vegetation on the gravity-driving term compared with the velocity model proposed by Huai et al. (2014). The upper inflection point z ₂ of the mixing layer is found above the top of the short vegetation canopy, while this point appears at the height of short vegetation in study of Huai et al. (2014). It is learned that the penetration length scale will change in different situations. The dense vegetation density (m ₁ = 778 plant·m⁻²) in this study may cause larger scale coherent vortices and thus penetrate to the tall vegetation layer, compared with the sparse vegetation density (m ₁ = 117 plant·m⁻²) in Huai et al. (2014); Figure 6 shows that the lateral vorticity reaches a peak at the height of short vegetation and decreases to zero near z = .075–.08 m, verifying the existence of the large-scale Kelvin-Helmholtz vortices. Moreover, according to Figures 4E, 8A in the study of Anjum and Tanaka (2020), the Reynolds stress turbulence model of double-layered vegetation flow indicates that the mixing region affected by the vertical coherent vortices contains the zone in short vegetation layer and tall vegetation layer, it is reasonable to assume that the upper inflection point of the mixing layer lies above the short vegetation canopy.

A slightly increasing tendency of the velocity in zone 1A (Figure 5), indicates the insignificant effect of the near-bed boundary friction. Recall from Eq. 2, the bed shear stress friction τ b e d 1 − λ 1 W d x was compared against the vegetation drag 1 2 ρ C D m 1 D U U D 1 2 W z 2 − δ w d x in zone 1A, which is given as: τ b e d 1 − λ 1 W d x 1 2 ρ C D m 1 D U U D 1 2 W z 1 d x = f b e d 1 − λ 1 4 C D m 1 D z 1 . (13)

For run YHS1, the λ 1 = 0.022 , C D = 1.68 , z 1 = 0.055 m , and the f b e d = 0.01 − 0.05 of the plexiglass plate was estimated by a Moody chart roughness height k s = 0.0025 m m (Zhang et al., 2020). Substituting the maximum value of the roughness parameter f b e d = .05 into Eq. 13 and yields a result of 2.80%. It means that compared with the vegetation drag the near-bed boundary friction is negligible in the momentum balance equation by comparing with the vegetation drag. Similarly, according to Eqs 2, 3, the sidewall friction term 2 z 1 τ w a l l d x and 2 H − z 2 τ w a l l d x were compared with the vegetation drag in zone 1A and zone 2A, respectively, as the following equations: 2 z 1 τ w a l l d x 1 2 ρ C D m 1 D U U D 1 2 W z 1 d x = f w a l l 2 C D m 1 D W , (14) 2 H − z 2 τ w a l l d x 1 2 ρ C D m 2 D U U D 2 2 W H − z 2 d x = f w a l l 2 C D m 2 D W . (15)

For run YHS1, f w a l l = 0.05 is selected to calculate Eqs 14, 15, and gets results of 1.06% and 3.19%, respectively. Even if the roughness conditions f _bed and f _wall being taken as the maximum value .05 in the sparsest vegetation arrangement run, the bed friction and sidewall friction reveal a small proportion relative to vegetation drag. Therefore, it is reasonable to assume that the boundary friction is negligible for all runs.

Eqs 6, 7 predict that the velocity is constant in zone 1A and zone 2A, but measured velocity has an increasing trend in zone 2A in Figure 7 and the model needs to be further improved. The velocity increasing is likely put down to the wake vortices generated by the cylinders, which are limited by the tall cylinders but not totally eliminated compared with the single-layered vegetation flow. In submerged vegetation flow, the wake region effect term can be approximately presented as a trigonometric function (Nikora et al., 2013): U W F = u * m 2 Π κ sin   2 π z 2 H . (16) where u * m is the shear velocity and is usually expressed as a formula of the maximum values of Reynolds stress which is usually found at the top of the vegetation layer (Järvelä, 2005; Nepf, 2012) i.e., u * m = − u ′ w ′ ¯ max , and the u * m applied here is calculated using the Reynolds stress at the short vegetation height, κ = .40 is the Von Kármán constant, П is Coles’s wake strength parameter and was suggested equal to .2 in uniform flows (Nezu and Rodi, 1986).

For the density of tall vegetation layer is sparser of Y arrangement than X, the increasing tendency of the velocity at zone 2A is more significant of runs YHS1 and YLS1, and the wake region effect is discussed of run YHS1. As shown in Figure 10, considering the wake term in the flow velocity model indeed can better reflect the experimental velocity characteristic in zone 2A. The Coles’s wake strength parameter changes from .2 to a value of .1, which means the resistance of the tall vegetation can significantly reduce the Coles’s wake strength parameter. The applicability of the wake parameter value to the velocity model of double-layered vegetation needs to be further explored.