Derivation of Canopy Resistance in Turbulent Flow from First-Order Closure Models

Quantification of roughness effects on free surface flows is unquestionably necessary when describing water and material transport within ecosystems. The conventional hydrodynamic resistance formula empirically shows that the Darcy–Weisbach friction factor f ~(r/hw) describes the energy loss of flowing water caused by small-scale roughness elements characterized by size r (<<hw), where hw is the water depth. When the roughness obstacle size becomes large (but <hw) as may be encountered in flow within canopies covering wetlands or river ecosystem, the f becomes far more complicated. The presence of a canopy introduces additional length scales above and beyond r/hw such as canopy height hv, arrangement density m, frontal element width D, and an adjustment length scale that varies with the canopy drag coefficient Cd. Linking those length scales to the friction factor f frames the scope of this work. By adopting a scaling analysis on the mean momentum equation and closing the turbulent stress with a first-order closure model, the mean velocity profile, its depth-integrated value defining the bulk velocity, as well as f can be determined. The work here showed that f varies with two dimensionless groups that depend on the canopy submergence depth and a canopy length scale. The relation between f and these two length scales was quantified using first-order closure models for a wide range of canopy and depth configurations that span much of the published experiments. Evaluation through experiments suggests that the proposed model can be imminently employed in eco-hydrology or eco-hydraulics when using the De Saint-Venant equations.


Introduction
The modeling of urban constructed wetlands requires routing an inflow hydrograph through vegetated canopies, where the vegetation may be emergent or submerged [1][2][3][4][5]. Such routing is commonly modeled by combining the continuity and the De Saint-Venant equations (SVEs) along the streamwise direction [6][7][8]. In this usual representation, the determination of the friction slope Sf necessitates a closure model, the subject of the work here.
can be used to link the bulk velocity for whole depth Ub to Sf, where n is Manning's roughness coefficient, and Rh is the hydraulic radius. To proceed further, a link between n and roughness element size r, given by n~r 1/6 , has been empirically proposed and is often referred to as the Strickler-scaling. This scaling describes several data sets reasonably as discussed elsewhere [11] but not over the entire range of ~r/hw [12], where hw is the water depth. Other versions replace Manning's formula with a Chézy equation where a coefficient C is used to represent roughness effects on the flow. However, these equations can all be made identical by linking their roughness parameterizations via   For steady and uniform flow at very high bulk Reynolds numbers, and upon further assuming the small-scale roughness (r/hw << 1) dominates the flow (left-side of Figure 1), traditional boundary layer theory leads to a mean velocity profile given as [13].
where κ = 0.41 is the von Kármán constant, z is the vertical distance from the channel bed, d is the zeroplane displacement, and z0 is the momentum roughness height that can be linked to n. Thus, the bulk velocity for flow over small-scale roughness (r/hw << 1) at high Reynolds number (i.e., for thin buffer and wake regions when compared to the log-region) can be linked to the log-law via [14].
where e is the Euler number. When hw >> d, η = 1/7-1/6 thereby recovering the aforementioned Strickler scaling for a certain range of large r/hw [12] for r~zo. The power-law approximation to the log-function in Equation (2) only applies for a certain range of hw/zo >> 1 discussed and delineated elsewhere [12,14]. However, when the roughness height is comparable to the flow depth, as may be encountered in aquatic vegetation covering a wetland or a river channel (right-side of Figure 1), the flow region can be decomposed into three layers according to the dominant vortical structures. Moving from the bottom to the flow surface vertically, the turbulent mixing-length is dominated by three vortical flow types, namely, von Kármán street (being spawned from the canopy element wakes), mixing layer (being formed by the fast flow above the vegetation and the slow flow inside the vegetation), and attached eddies far from the vegetated surface as expected from canonical roughwall boundary layer theory [15]. Support for this 3-layer representation is as follows: (1) In the deeper vegetation layer (z/hv < 1), spectral analysis reveals that the energetic motion is dominated by von Kármán vortex streets with size Lv, which is proportional to the vegetation stem diameter D, where the mixing length is strictly dependent on the vegetation element diameter and independent of the local velocity and canopy density; (2) well above the canopy (z/hv >> 1), the flow resembles a canonical rough-wall boundary layer, where the characteristic vortex size LCBL scales with the distance from the boundary (or zero-plane displacement) as expected for attached eddies; (3) for the zone near the canopy top (z/hv ≈ 1), the flow is dominated by mixing-layer eddies [16] but occasionally is disturbed by attached eddies, although these types of eddies do not coexist in space. Simulations and experiments provide substantial support for this hypothesis as reviewed elsewhere [17]. The dependence of the mixing length on the typical vortex size for each layer is summarized as [15].
for von Kármán streets (5) where l is a canonical mixing length, and L is the vortex length scale. Figure 1 makes clear why the flow structure within the vegetation is considerably more complicated than small-scale roughness. The velocity profile also changes considerably from its canonical log-law (or law of the wall). Many approaches can be used to model u(z), such as the numerical models using first-order or higher-order closure schemes [14,18], and analytical models [19,20] with different expressions of Reynolds stress closures. Moreover, the velocity profile for the entire flow depth is complicated and varies with the vortex development in different layers [21][22][23][24]. These 'distortions' to the classic logarithmic mean velocity profile and the possibility of deriving a general expression for friction to be used in operational models (such as the SVE) for submerged vegetation frame the scope of this study. The main finding is the dependence of the friction factor on two dimensionless groups related to canopy submergence and canopy length scale. The dimensionless groups are determined using first-order closure model calculations that utilizes the mixing length scales in Equations (3)-(5). The bulk velocity as well as the friction factor derived using this approach are compared to a large corpus of published experiments where the parameters needed for experimentally computing the two dimensionless groups can be inferred. The usage of the first-order closure model to infer the two dimensionless groups is advantageous because it allows explorations of parameter combinations and flow conditions not covered by current experiments.

Vegetation Resistance
In this section, the vegetation resistance per unit bed ground is given first, which requires two flow-rate parameters describing the effects of vegetation elements on the flow: The drag coefficient Cd (discussed in Section 2.1.1) and the choice of the reference velocity within the vegetated region of the flow Uv (discussed in Section 2.1.2). Then the friction factor of vegetation is given in Section 2.1.3 based on the approach of Darcy-Weisbach friction using the results of Sections 2.1.1 and 2.1.2. The novelty of the proposed approach here lies in their joint specification via two dimensionless groups that accommodate the vortical structures featured in Figure 1.
For per unit bed area, the resistance by vegetation in submerged conditions is dominated by form drag, which can be calculated as [21,[25][26][27] where Cd is the overall vegetation drag coefficient, which is a function of the flow state quantified by a Reynolds number (to be discussed later); m is the vegetation density determined by the number of stems per unit bed area; D is the frontal width of the canopy; and Uv is a reference velocity representing the flow in the vegetation section. There are several plausible velocities to represent Uv, and those are discussed later. Beyond the vegetation properties (m, D, hv), Equation (6) requires two flow-rate parameters describing the effects of vegetation elements on the flow: The drag coefficient Cd and reference velocity of vegetated region Uv.

Drag Coefficient
The drag coefficient Cd, which codifies the interaction between the flow state and the morphology and density of the canopy, has been widely studied [28,29]. When neglecting the interaction between adjacent wake elements of a vegetation array, the standard expression [7,9,30] for the drag coefficient of an isolated cylinder-type canopy for the range of Red between 0.02 and 2 × 10 5 is: where a local element-related Reynolds number can be expressed as: For very large Reynolds number (>2 × 10 5 ), Equations (7-9) predict a constant Cd,iso ≈ 1.2. In most cases, the interaction between wakes of elements cannot be overlooked; thus, the effect of the canopy array is not the same as that of an isolated one [9,31]. For this reason, the drag coefficient of a vegetation array has attracted significant research interest [14,30,[32][33][34][35][36].
Some studies focused on variations of the cross-sectional averaged drag coefficient along the streamwise direction. For steady non-uniform flow through a uniform-arrangement of emergent vegetation, the cross-sectional averaged Cd exhibits a parabolic-shape along with Reynolds number [9], indicating that the drag effect initially increases and then decreases along with the Reynolds number. When adding rainfall on nonuniform flow in laboratory experiments [7], the rainfall also changes the features of the cross-sectional drag coefficient, which changes the parabolic shape of the Cd-Re relation into a monotonous one where the drag coefficient decreases along with increasing Reynolds number during large rainfall events.
Others have investigated the variation of the drag coefficient along the vertical direction. Poggi et al. [15] experimentally considered the effect of vegetation density on the canopy sub-layer turbulence. Vertical variation of the drag coefficient within a rod canopy was obtained according to a simplified mean momentum balance equation for the streamwise direction. A linear empirical function between Cd-Re was proposed on the basis of measured profiles of Reynolds stress and the velocity from laboratory experiments. Naturally, this linear relation applies for a restricted range of water depth and Red.
For the volume-averaged drag coefficient, Tanino and Nepf [37] conducted experiments in two Plexiglas recirculating flumes with randomly arranged cylindrical maple dowels. Their experiments also resulted in a linear function of the volume-averaged Cd-Re. The monotonic decline trend was also found in other experiments. Cheng and Nguyen [38] assembled numerous experimental results [23,37,[39][40][41][42][43] and reported a monotonic decline in drag coefficient with an increasing vegetationrelated Reynolds number Rev. An expression that fits their data and incorporates the fractional volume concentration is given by: where the vegetation-related local Reynolds number is: where ϕ is area concentration of vegetation stems. Again, for very large Red (>2 × 10 5 ) and small ϕ, Rev reduces to (π/4) Red and Cd,array ≈ 0.8 (a constant).

Reference Velocity
Because the drag coefficient only senses the velocity within the vegetation, it becomes necessary to separate this flow from the bulk flow in models for f. Several plausible velocities are now introduced to represent this velocity.
(1) When considering the maximum flow width of channel B, the bulk velocity in the vegetation layer can be calculated as [39,44,45].
where Qv is flow rate for vegetation layer.
(2) In view of the blockage of the vegetation elements, the pore velocity in the vegetation layer can be the reference velocity, which is calculated based on a spatially averaged value as [7,9,37,38,43,46].
where the effective width Be=B(1-ϕ) is used instead of the maximum flow width B.
(3) In a classical staggered array, the constricted velocity is: with a characteristic constricted width Bc=B(1-D/Ls,stag), and Ls,stag is the spacing distance defined by Etminan et al. [47]. The work of Etminan et al. [47] also showed that for very large Red (>2 × 10 5 ) and small , Cd ≈ 1 (a constant) for a staggered cylinder arrangement.
(4) In view of the flow structure in turbulence, the separation velocity can also serve as a characteristic velocity, which is expressed as [47].
where kp is a kinetic energy of the sub-grid scale for a Smagorinsky eddy-viscosity model.
The spatially averaged pore velocity is adopted here (i.e., Uv = Uv,pore) due to the diverse vegetation arrangements covered by the experiments analyzed here.

A Friction Formula
Adopting a Darcy-Weisbach friction factor equivalence by imagining the vegetation layer as producing the same f as a thin roughness layer, the shear stress dominated by vegetation is expressed as Equation (6), which can be used to link vegetation drag per unit bed area to f as [48,49].
where the subscript v denotes the mean friction factor induced by vegetation; and Lc is the adjustment length scale, which is expressed as: The ratio of velocity in the vegetation layer to that for the bulk flow depth is a key component in estimating fv. Define Us as the depth-averaged velocity for the free water layer overlying the vegetation, and ∆U=Us-Uv is the bulk velocity difference between the surface and vegetation layers, then the depth-averaged velocity for the whole depth Ub can be expressed as: It directly follows that Uv/Ub in the friction factor of vegetation expression (Equation (17)) can be made as a function of submergence hv/hw and ∆U/Uv via where ∆U/Uv indicates the blocking effects of the vegetation on the flow. Previous results [49] showed that fv can be linked to two dimensionless groups, namely, submergence and drag length scale, which are defined as Combining Equations (20)- (22), the fv (Equation (17)) becomes a function of three terms: α, β and ∆U/Uv. A formula linking ∆U/Uv to α, β can be derived using a first-order closure model thereby allowing the friction factor to be uniquely varying as a function of α and β, i.e., fv = fv (α, β).

First-Order Closure Models
The streamwise momentum equation based on the double-averaging method yields [21,24].
where Du/Dt is the total or material derivative, ρ is the water density, ν is the kinematic viscosity, So = sinθ is the bed slope, g is the gravitational acceleration, z is the height above the bed. Term u w z     is the spatially averaged Reynolds stress whereas the term u w z     is the dispersive stress arising from spatial correlations in the time-averaged velocity field. It has been shown by Poggi et al. [50] that the dispersive stress is less than 10% of the Reynolds stress for mDhv > 0.  [15,49].
The drag force resulting from the presence of a canopy can be approximated by the quadratic law and is given as: For the Reynolds stress, a conventional K-theory closure yields.
where Km is the eddy diffusivity for momentum and is impacted by the vortical structure dominating the various layers as shown in Figure 1. By substituting Equations (25) and (26) into (24), the momentum equation becomes: Two boundary conditions, a model for diffusivity Km and a model for the drag coefficient Cd are required to solve the above second-order ordinary differential equation.
The eddy diffusivity can be modeled using mixing-length arguments as: where leff is the effective mixing length parameterized as a function of the established vortex sizes of the different layers described in Figure 1.
In the canonical boundary layer, the classical attached eddy hypothesis provides an estimate for: where d is the zero-plane displacement height. The zero-plane displacement may be evaluated in multiple ways but it is commonly determined from the center-of-pressure method [51] given as: The vegetation layer consists partly of mixing layer eddies and von Kármán vortex streets within the entire canopy volume but most effective in the bottom layers of the canopy. For the Kármán vortex street layer, the mixing length and the size of the characteristic vorticity are approximately equivalent. Hence, for the entire vegetated layer, it is assumed that: where  is a constant determined by enforcing a continuity (but not smoothness) of the effective mixing length across the entire flow depth. Enforcing a continuity in mixing length yields.
Selecting the drag coefficient is also significant to first-order closure modeling and the formulation based on the array configuration is used.
The data sets span wide-ranging flow and canopy properties. The effective width of the canopy layer was determined as B (1-ϕ), which considers the volume occupied by vegetation [9,37]. The entire range of Reynolds number (Red=UbD/ν) for these data ranged from 61 to 9936. The Froude number (Fr=Ub/(ghw) 1/2 ) ranged from 0.0045 to 0.5649 (sub-critical for all runs). The vegetation drag coefficient (Cd, array of Equation (11)) ranged from 0.84 to 6.35 with the averaged value 1.28 for all the data points here. A summary of the basic information for the experiments analyzed here is featured in Table 1.  Figure 2 compares the bulk velocity modeled using the first-order closure model (solving momentum Equation (27) along with Equations (28) to (32)) and the measured bulk velocity. The drag coefficient Cd,array (calculated by Equation (11)) was adopted in this model, and the value is shown in Table 1. The results show that using Cd,array derived from numerous experiments reasonably predicts measured bulk velocity when employed with the effective mixing length featured in Figure  1. The modeled bulk velocity was determined from the first-order closure by using the computed mean velocity profile via: This agreement lends confidence to using the first-order closure model to assist in the determination of the two aforementioned dimensionless groups for term ∆U/Uv.

First-Order Closure Model Runs
In this section, term ∆U/Uv reflecting the blocking effects of vegetation on the flow is investigated. A link between ∆U/Uv and the two dimensionless groups α and β is to be established using the first-order closure model results. In Section 4.1, a scaling analysis is first conducted between ∆U/Uv and α while fixing β. In Section 4.2, a similar scaling analysis is to conducted between ∆U/Uv and β while fixing α. Finally, In Section 4.3, a general relation between ∆U/Uv and α , β is proposed by combining results from Sections 4.1 and 4.2. Last, an empirical formula for ∆U/Uv to be used in the determination of f is proposed and comparisons with bulk velocity and friction factor measurements are conducted.

Scale Analysis with Submergence
When fixing β, ∆U/Uv must vary only with α. As customary with scaling analysis, it is assumed that: When the height of surface layer hs = hw − hv = 0, indicating an emergent condition, a ∆U=0 is recovered.
Here, the first-order closure model is adopted to explore the coefficients in Equation (34). We only focus on submergence and keep other parameters fixed.  Figure 3 shows the best fit curve to results obtained from the first-order closure model.

Scale Analysis with Vegetation Attributes
The blockage effect term ∆U/Uv also increases with vegetation drag coefficient originating from the frontal area per unit bed area. Thus, an empirical formula can be given as follows based on a scaling analysis.
When no vegetation exists (that is, m, D, or hv equals zero), ∆U=0 is recovered. Equation

Expression for the Combined Influences of Submergence and Vegetation
A general relation between ∆U/Uv and α, β is now proposed by combining Equations (34) and (35) to yield.
To cover the entire range of α (submergence attributes) and β (vegetation attributes), the basic conditions (submergence and vegetation attributes) of the experiments were used in the first-order closure model. The results from the first-order closure model are shown in Figure 5. It can be seen that Equation (36) reasonably collapses the first-order closure model results. The model in Equation (36) and its use in the scaling analysis leading to Equation (17) constitutes one of the main novelties of the work here. Specifically, the friction factor formula can be obtained as: with best fitting parameters c1 = 1.7237, c2 = 0.8545, and c3 = 0.4944. Hence, it follows that: where the depth-averaged pore velocity can be obtained from the mean momentum balance to yield [48,62].

Model Validation
In this section, the derived expression for ∆U/Uv is used to model the friction factor fv and bulk velocity Ub and then compare against laboratory experiments described in Table 1. The comparisons between model calculations (Equations (37)- (39)) and data are featured in Figures 6-8. It is to be noted that α, β here are determined from the first-order closure model results and thus c1 = 1.8629, c2 = 0.7909, and c3 = 0.5137 are simply a summary of those model calculations. No data were used in the determination of α, β or coefficients c1 to c3. Hence, the comparisons with laboratory experiments in Figures 6-8 can be treated as validation.
Comparisons between the calculated (using Equation (36)) and measured ∆U/Uv are presented in Figure 6. Almost all the data points collapsed on a 1-to-1 line except for six points when ∆U/Uv was small (submergence was minimal), and they were affected by many other factors. The flow in the surface layer may not have been fully developed in those cases. From Figure 6, agreement between modeled and measured is acceptable for ∆U/Uv > 0.3.

U/U v by 1st order closure model
Although comparisons between the calculated and measured fv show scatter, the modeled bulk velocity reasonably predicts the measurements. The Mean Squared Error (MSE) is used here to evaluate deviations between modelled and measured results and is given as:  (40) where Nsample is the number of the data sample. For ∆U/Uv, the comparison between modelled and measured variables yields MSE = 0.0894, and maximal departure from modelled to measured is 1.5565.
From another point of view, it is to be noted that Figure 5 presents a comparison between modelled (using Equation (36)) and results from first-order closure model calculations, whereas Figures 6 shows the comparison between modelled (also using Equation (36)) and measurements summarized in Table 1. Similarity of these two figures also implies that the first-order closure model also recovers the experiments reasonably. For fv, the comparison between modelled and measurements yields MSE = 0.2706 and maximal departure of 2.8152. For bulk velocity Ub, comparison between model and measurements yields MSE = 0.0041 m 2 /s 2 , with maximal departure between model and measured being 0.3581 m/s.

U/U v measured
In summary, using a first-order closure model whose results are summarized in Equation (36), the "standard velocity difference" between the free water layer and the vegetation layer ∆U/Uv can be modeled when ∆U/Uv > 0.3. Using Equations (37) and (38) that are derived from Equation (36), the f v,measured U b,measured friction factor and bulk velocity can be reliably determined and used for modeling overland flow via the SVEs.

Conclusions
In comparison with small-scale roughness values, canopies introduce additional complications and length scales above and beyond r/hw such as canopy height, arrangement density, frontal element width, and drag coefficient. To link those length scales to the friction factor, scaling analysis aided by first-order closure model calculations were adopted. It was shown that the Darcy-Weisbach friction factor for canopies varies with two dimensionless groups, namely, canopy submergence and a dimensionless canopy length scale. These two dimensionless groups were then determined from a first-order closure model calculation that explicitly considered vortical sizes within and above the various canopy zones. Comparison between experiments and calculated bulk velocity and friction factor suggested that the proposed link between Sf and bulk velocity employed in SVEs can be used for describing flow through vegetation.