Hybrid modeling of flows over submerged prismatic vegetation with different areal densities

In the modeling of vegetated flows an efficient approach is to solve the Double Averaged Navier Stokes (DANS) equations, which are obtained from the spatial and temporal averaging of the Navier Stokes equations. The resistance effects of vegetation are modeled by a drag force density term with an empirical bulk drag coefficient, as well as by turbulence terms characterized by a length scale parameter. These empirical parameters are dependent on the areal density of vegetation and the spatial distribution pattern of vegetation elements and have seldom been studied. In this work the effect of spatial distribution of vegetation elements on the bulk drag coefficient is investigated by computing explicitly the flows around the vegetation elements. The results show that the bulk drag coefficient increases with the longitudinal vegetation element spacing, decreases with the lateral vegetation element spacing and can have multiple values for a given areal density of vegetation. In the DANS model the vegetation induced turbulence is simulated by a novel k- (cid:2) type model embrac-ing an empirical length scale parameter. The length scale parameters are calibrated against previous experiments and a new set of experiments with high areal density of vegetation. The DANS model is subsequently verified by two cases with the same vegetation density and different distribution patterns.


Introduction
The growth of vegetation in open channels contributes to the maintenance of ecological balance there as well as the sustainability of the environment. Vegetation will increase turbulence of flow and enhance mixing of solute and suspended substance. It however also increases hydraulic resistance and reduces channel flow carrying capacity. It is important to determine the flow resistance generated by vegetation for the engineering design of vegetated channels and the flow characteristics for the assessment of the environmental impact.
Hydraulic resistance in a vegetated channel is due to the drag force predominately generated by the pressure variation on the surfaces of vegetation stems. For generality the pressure drag is represented by a dimensionless drag coefficient multiplying the dynamic pressure. Within a vegetation zone the drag of a stem will be affected by the neighboring stems due to wake interference, and the bulk drag coefficient will also be affected. Tsujimoto, Kitamura, and Okada (1991) conducted measurements and modeling of flow and turbulences generated by simulated vegetation, and reported good agreement between experiment measurements and CONTACT Chi-wai Li cecwli@polyu.edu.hk the numerical results using a modified k-with dragrelated source terms. Nepf (1999) studied the conversion of mean kinetic energy into turbulent kinetic energy by simulated vegetation and developed a model to describe the drag and turbulence of flow through emergent vegetation. Lopez and Garcia (2001) carried out spatial and temporal averaging of the conservation laws to develop two-equation turbulence models and used them to analyze the mean flow and turbulence structures of open channel flow through submerged vegetation. Stone and Shen (2002) carried out extensive laboratory tests to study the hydraulic resistance of flow in channels with submerged or emergent cylindrical stems, and developed formulas relating hydraulic roughness to flow depth, areal density of stems, stem length and stem diameter. For flows over submerged vegetation, within the vegetation zone the flow is governed by the drag resistance induced by the vegetation stems. Above the vegetation the flow is governed by the shear resistance at the interface of the vegetation zone and clear water zone (e.g. Baptist et al., 2007;Ghisalberti, 2009). In the modeling of submerged vegetated flows the double-averaged Navier Stokes (DANS) equations (e.g. Li & Li, 2015;Lopez & Garcia, 2001) or the Reynolds Averaged Navier Stokes (RANS) equations (e.g. Tsujimoto et al., 1991;Zhang, Li, & Shen, 2010) are commonly used instead of directly computing the flows around each piece of vegetation. The double-averaging procedure generates some terms that needed to be parametrized. The parametrization process generates empirical coefficients, including the bulk drag coefficient of vegetation and a parameter accounting for a shift up of the velocity profile (similar to the zero-plane displacement parameter, e.g. Baptist et al., 2007).
The dependence of bulk drag on vegetation density for cylindrical stems has been shown experimentally by Stone and Shen (2002) for staggered arrangement of stems, and by Tanino and Nepf (2008) for random arrangement of stems. For blade-type vegetation with rectilinear arrangement of vegetation elements, Busari and Li (2016) showed that the bulk drag is not only dependent on vegetation density, but also on the longitudinal and lateral separation of adjacent stems. The generalization of the results by using an empirical formula is not straightforward as there are many parameters involved, such as the cross-sectional shape of the vegetation element and the relative location of adjacent elements. Instead, a numerical modeling approach can be used to compute the flow around individual vegetation elements subjected to wake interference. Koch and Ladd (1997) computed flows through cylinder arrays for Reynolds number up to 180 using the lattice-Boltzmann method. Stoesser, Kim, and Diplas (2010) carried out 3D Large Eddy Simulation (LES) of flows through staggered arrays of emergent cylinders with three values of areal density of vegetation. The turbulence structures were studied in detail. The pressure drag and skin friction drag were found to contribute to over 98% of the total energy loss for vegetated flows with moderate to high areal density. The results indicated that the three-dimensional effects on drag characteristics may be very small.
For submerged vegetated flows, the bulk drag of vegetation restricts the amount of flow passing through the vegetation. The shear at the vegetation top produces coherent motion and hydraulic resistance (Ghisalberti & Nepf, 2002;Nezu & Sanjou, 2008). An inflection occurs on the profile near the vegetation top (Ghisalberti, 2009). Li and Li (2015) developed a DANS model with a modified SA turbulence closure. The turbulence length scale is shown to be dependent on the areal density of vegetation. Lopez and Garcia (2001) developed a k-model with drag related sink terms to account for vegetation generated resistance and turbulence. In the model the effective turbulence dissipation time scale is assumed equal to the intrinsic turbulence time scale. Uittenbogaard (2003) developed a k-model with a modified time scale term to account for vegetation generated turbulence. Derived from the time scale terms, a length scale accounting for the penetration of shear turbulence from the clear water zone at vegetation top into the vegetation zone was proposed. Souliotis and Prinos (2011) used the VARANS (volume averaging of the RANS equations) approach and compared four different turbulence models in the simulation of open channel flows with submerged vegetation.
In the modeling of vegetated flows using the DANS approach the bulk drag coefficient and the shear induced turbulence length scale parameters are the major sources of modeling errors. In this work the effects of areal density of vegetation and the distribution pattern of vegetation elements on these parameters are investigated numerically and supported by experiments. Firstly the bulk drag coefficient is determined separately using a 2D RANS model with a high resolution grid such that the flow around individual vegetation elements is computed and the velocity and pressure fields are obtained. The wake interference effects among stems are studied through simulation of cases of vegetation with blade-type cross-section or cylindrical stems arranged in regular or irregular array patterns, and with the areal density of vegetation elements varied by changing the longitudinal spacing or lateral spacing of vegetation elements. Secondly a DANS model with a k-epsilon type turbulence closure is developed with the additional source terms related to the areal density of vegetation and bulk drag coefficient. The computed bulk drag coefficient from 2D RANS model is input into the DANS model to predict vegetated flows.

Governing equations
To derive the DANS equations, the RANS equations are spatially averaged over a domain of the order of the vegetation stem spacing much larger than the vegetation generated turbulence length scales and much smaller than the characteristic large scale of the flow. For constant porosity, the DANS equations can be written as follows (Nikora et al., 2007).
where x i ( = x, y, z) are the coordinates in longitudinal, transverse, and vertical directions respectively; u i (= u , v , w ) are the double-averaged velocity components in x, y and z directions respectively; t = time; ρ = density of fluid; p = double-averaged pressure; F i (= F x , F y , F z ) are the double-averaged resistance force components per unit volume induced by vegetation in x, y and z directions respectively; g i ( = g x , g y , −9.81 m/s 2 ) are the components of the gravitational acceleration; ν m = molecular viscosity; − u i u j = double-averaged correlation term of temporal fluctuation of velocity = space-averaged Reynolds stress; and − ũ iũj = space-averaged correlation term of spatial fluctuation of velocity = form-induced stress or dispersive stress.
The resistance force due to vegetation arises from the surface integral of the normal pressure and lateral shear stress on the vegetation stems. Generally the surface integral of the shear stress is small. The surface integral of the pressure and shear on vegetation can be determined by a quadratic friction relationship. The spatially average force per unit volume within the vegetation domain is obtained by where C d = bulk drag coefficient; b = stem width; N = areal density of vegetation (number of stems per unit area, in 1/m 2 ). The determination of the bulk drag coefficient will be discussed in the next section. The space-averaged Reynolds stress can be parameterized in the following way as proposed by Pedras and de Lemos (2001).
where ν t = eddy viscosity; δ ij = Dirac delta function; k = turbulence energy. The last terms of the equations can be absorbed into the pressure term. Some extended k-models have been proposed for vegetated flows. Examples include Baptist et al. (2007), Lopez and Garcia (2001), and Defina and Bixio (2005). In the present work the following form is proposed: The last term of the k equation describes the turbulence energy generation due to vegetation. The last term of the equation describes the turbulence dissipation rate due to vegetation, where L p = length scale of turbulence energy dissipation. The set of standard constants is C μ = 0.09, C 1 = 1.44, C 2 = 1.92, σ k = 1.0, σ ε = 1.3 (Rodi, 1984). The distinctive feature of the model is that L p is related to the vegetation density and drag coefficient. The determination of L p is discussed in section 4 and the empirical expression of L p is given by Equation (6).
The numerical model developed is based on Lin and Li (2002) and Chen, Li, and Zhang (2008). The computational domain with variable free-surface and irregular bottom topography is transformed to a rectangular prism by σ -coordinate transformation. The governing equations are solved on a non-uniform rectilinear grid by a split operator finite difference method. Proper boundary conditions are imposed at the inlet, outlet, surface, bottom, and lateral boundaries.

Determination of drag coefficient
The drag coefficient of an isolated cylinder over a wide range of Reynolds number has been determined (e.g. Munson, Young, & Okiishi, 2002). For vegetation with rigid stems it is often simulated by a group of cylinders or rectangular plates. The bulk drag coefficient will depend on the shape of individual vegetation elements and the areal density of the stems, as well as the distribution pattern of the vegetation elements (Busari & Li, 2016;Stone & Shen, 2002;Tanino & Nepf, 2008).
As the wake interference effect among vegetation elements are complicated and there are many variables (shape, areal density, lateral and longitudinal spacing between adjacent vegetation elements) needing to be considered, a general expression of the bulk drag coefficient for arrays of vegetation elements is difficult to obtain. The engineering approach is to assign a typical value for the bulk drag coefficient, usually the value for an isolated vegetation element, assuming the wake interference effect is not important. For an array of cylinders, a value of C d ∼ 1.2 is widely used. For blade-type vegetation element, a value of C d ∼ 2 is used. These values however are not accurate for vegetation with high areal density. Tanino and Nepf (2008) found that C d > 4 for an array of cylinders with solid area ratio greater than 0.3. Busari and Li (2016) found that C d > 6 for an array of blade-type elements with narrow lateral spacing between adjacent elements.
In the present work a RANS numerical model is used to compute the flows around vegetation elements. The pressure field and velocity field are determined and the drag coefficient is directly calculated. The vegetation elements are emergent and the wake interference behavior is assumed two-dimensional. The three-dimensional effect is considered not significant. The computation is very efficient and many cases can be simulated in a short time.
In the first set of numerical simulations the vegetation element is of blade-type. The elements are arranged with a rectangular grid pattern. For each run the velocity of flow, lateral and longitudinal spacing of vegetation pieces are varied. As the number of runs is large, the property of symmetry is utilized. The computational domain covers a region enclosing nine pieces of vegetation. The flow is periodic at the inlet and outlet, as well as at the right and left sides. The vegetation is emergent and the flow is assumed two-dimensional. The schematic diagram of flow over vegetation is shown in Figure 1. Figure 2 shows the grid layout and imposed boundary conditions of the computational domain.
A body force is introduced into the domain such that the areal average velocity is equal to the specified value. The drag force on each piece of vegetation is computed by summing all the pressure forces and shear forces acting on the surface of the piece of vegetation: a turbulence model is required. Two turbulence models have been tried: the standard k-model and the k-ω SST model. Grid convergence tests have been carried out by using three different grids. The final grid chosen is 210 × 210 shown in Figure 1. The grid size is variable and is finer around the pieces of vegetation and coarser at the middle region of two adjacent pieces of vegetation. The model is run until dynamic steady state is reached.

Shielding effects
The flow induced drag force on an object will be hindered by the presence of another object upstream. This is referred as the shielding effect.  (Busari & Li, 2016), as shown in Figure 3. For the use of the k-model the difference in the computed and measured bulk drag coefficients are within 10%. For the use of the k-ω SST model the difference between the computed and measured values is slightly larger. Further confirmation of the decreasing trend C d with the increasing of S y is done by conducting more cases of computation with values of S y changes to 0.0125 and 0.025 m. Details of the parameters are shown in Table 1.
The vortex shedding of flow around a piece of vegetation is maintained for the case with the ratio S x /S y  being not too large (< 2), as shown in Figure 4(a). At a higher ratio (S x /S y > 2) the strength of the vortex shedding is reduced (Figure 4(b)). At smaller spacing S x each piece of vegetation is affected by the wake generated by the upstream piece of vegetation. The pressure on the upstream face thus is lower and the resulting drag force is reduced. The large value of the bulk drag coefficient at large ratio of Sx/Sy is mainly due to the narrow lateral spacing between adjacent vegetation elements, which induces large energy loss.

Channeling effect
The channeling effect refers to the increase in drag on an individual piece of vegetation when the transverse spacing of adjacent pieces of vegetation decreases. To illustrate this, three cases of difference transverse spacing and fixed longitudinal spacing are simulated. The longitudinal spacing S x is fixed at 0.02 m and the lateral spacing used is 0.02, 0.05, and 0.08 m. The computed results show that the drag coefficient increases with the decrease of lateral spacing. The difference between the computed results and measured data (Busari & Li, 2016) is also around 10% ( Figure 5). Furthermore, additional tests with different value of Sx changed to 0.01 and 0.03 m are conducted to support the above observations. Details of the parameters are shown in Table 1. For a given discharge rate or areal mean velocity of flow, the decrease of the transverse spacing of adjacent vegetation elements increases the velocity in between two adjacent pieces of vegetation. The pressure exerted by the flow on the upstream surface of a piece of vegetation is increased and the resulting drag on the piece of vegetation is increased. As the reference velocity used is the average velocity the drag coefficient is increased due to the increasing drag. Vortex shedding occurs for both cases as shown in Figure 6. The strength of the vortex shedding is stronger for the case with narrower spacing.

Effects of shape and distribution pattern
For many species of vegetation the cross-sectional shape of a stem is close to a circle. The effects of circular shape and the layout of stems are studied numerically in this section. The circular stems are arranged in a rectilinear  pattern or in an irregular pattern (Figure 7). The computational domain covers a region enclosing nine cylinders. The grid used is sufficiently fine with 271 × 271 grid points. The Reynolds number is 1000 with the velocity taken to be the average velocity over the pore space between stems. The solidity ratio of the vegetation ranges from 0.03 to 0.35.
For the stems of irregular distribution pattern, the drag coefficients for individual stems are different. The bulk drag coefficient is computed by averaging all the stem drag coefficients. The results are shown in Figure  8. It can be seen that the bulk drag coefficient increases with the areal density of vegetation and is in agreement with experimental results on random arrays of cylinders (Tanino & Nepf, 2008) and staggered arrangement of cylinders (Stone & Shen, 2002). The agreement is very good for cases with lower areal density of vegetation. The discrepancy is larger for the case with the highest areal density of vegetation. This can be explained by the fact that for high areal density of vegetation, the flow encounters the next downstream cylinder in a shorter distance, and will be affected more significantly by the exact location of the downstream cylinder.
For the stems with rectilinear distribution pattern, the shielding effect is significant for cases with high areal density of vegetation. The drag coefficients for individual stems are the same, since every stem is subjected to the same wake interference effects. As shown in Figure 8, the bulk drag coefficient also increases with vegetation density. The magnitude of the bulk drag coefficient is smaller than that from the corresponding experiments on stems of random or staggered patterns. These results show that the distribution pattern can have significant effect on the bulk drag coefficient.
The computed streamlines show that there is a clear path for the case with rectilinear array of cylinders ( Figure 9). This reduces the resistance to the flow and results in a lower bulk drag coefficient. The effect is apparent for cases with high areal density of vegetation as the flow will be obstructed by the downstream cylinder in a short distance.
Comparing Figures 8 and 3, it can be observed that the bulk drag coefficient for arrays of rectangular plates is generally larger than that for arrays of cylinders. This is mainly due to fact that the drag coefficient of an isolated rectangular plate ( ∼ 2) is larger than that for an isolated cylinder ( ∼ 1.2).
The results on the dependence of drag coefficient on the spatial distribution pattern of vegetation are not exhaustive. They can be used for cases with similar spatial distribution pattern and flow conditions. The methodology of using a CFD (Computational Fluid Dynamics) approach to compute drag coefficient however is general and has been shown to work properly in the present study. For cases not covered in the present study, the CFD approach can be used to determine the drag coefficient.

Calibration of turbulence length scale by experiments
At the vegetation top a shear layer is formed due to the retarded flow inside the vegetation. Coherent vortices will occur and the vertical exchange of momentum causes the large eddies at the vegetation top to squeeze into smaller scale eddies within the vegetation region.
The turbulence length scale L p should be related to the geometrical properties of vegetation, including the vegetation height h v , the stem width b, the stem spacing ∼ √ N, and the shape and distribution pattern of the stem which is reflected in the value of C d . L p can be interpreted as a kind of turbulence eddy intrusion depth. It will approach the vegetation height if the vegetation density is very low, and is expected to approach zero if the vegetation density is very high. The following empirical expression is proposed: where f v = C d Nbh v , γ and δ are empirical parameters. The shape of the curve is similar to that of Raupach (1994) in the determination of the displacement height. The appropriateness of the functional form and the value of the empirical parameters are then determined. There are some available experimental data suitable for model calibration, including Dunn, Lopez, and Garcia (1996), Nezu and Sanjou (2008), and Zeng and Li (2014). The list of relevant parameters for the cases is shown in Table 2. The range of f v is from 0 to 6.  To provide data for vegetation with high value of f v a new set of experiments was conducted. The set of experiments is an extension of the experiments on emergent vegetated flows by Busari and Li (2016), with the water depth increased such that the vegetation becomes submerged. The experiments were carried out in a tiling flume of width 0.31 m, depth 0.40 m, and length 12.50 m. The flow is generated by a pump connecting to a recirculation system. A flow straightener using honeycombs was installed at the entrance to streamline the flow. An electromagnetic flowmeter was installed in the flow return pipe to record flow rate. The longitudinal water surface profiles were measured by Vernier point gauges with ±0.5 mm accuracy. The vertical profiles of mean velocity and Reynolds stress at the centerline location 0.3 m upstream of the simulated vegetation domain were measured using a 3D side-looking ADV (Acoustic Doppler Anemometer).
The vegetation patches are of width 0.3 m, with length varying from 0.75 to 2.4 m; the higher the areal density of vegetation, the shorter the length will be. Arrays of semi-rigid cable tile blades are used to simulate vegetation. The cable tile blades are of height 0.25 m, width 0.0075 m, thickness 0.0017 m, and are fixed on a PVC board in the flume with a regular rectilinear grid pattern ( Figure 10). The blades are semi-rigid and their deflection is small under the laboratory flow conditions. The value of f v ranges from 2 to 36.
For each selected case the DANS model is run with the turbulence length scale optimized such that the computed velocity profile is closest to the measured profile. The driving force (body force) per unit mass is set equal to the measured surface water surface slope. The results of the best-fit values of L p /h v against the densityvegetation-density parameter f v is shown in Figure 11. The proposed equation (Equation [6]) is also plotted on Figure 11. It can be seen that the equation fits the data fairly well with δ = 0.3, γ = 0.4. The results show that the ratio L p /h v decreases rapidly from 1 to a value of 0.4 for the drag-vegetation-density parameter ranging from 0 to 1. It then decreases slowly towards 0 for f v approaches infinity. At this range of drag-vegetation-density parameter the vegetation retards the flow significantly and displaces most of the flow upwards, causing a shift-up of the velocity profile towards the vegetation top. It should be noted that the values of the empirical parameters in  To determine the drag coefficient the 2D model using OpenFoam is utilized. The arrangement of the vegetation elements for the two cases is shown in Figure 12. The snapshots of pressure field and streamlines are also shown in Figure 12. For the case with narrower lateral spacing, the flow resembles an orifice flow with significant flow contraction and large head loss, resulting in a large drag coefficient of 10.0 and the value of f v = 15.0. The phenomenon of vortex shedding is suppressed. For the case with wider lateral spacing, the flow around a vegetation element is less interfered by the surrounding elements. Vortex shedding still occurs and the resulting drag is smaller with an average value of 3.1 and the value of f v = 4.7. For the turbulence length scale, both cases fall into the high range of areal density of vegetation. The value of the parameter L p /h v is 0.25 for case with lower of C d , which is higher than the value of 0.17 for the case with higher C d .
In the simulation, the effect of surface slope is represented by the body force per unit mass. The water depth is fixed at the measured value. At steady state the body force is balanced by the vegetation generated resistance force. The computed vertical profile of the longitudinal velocity and Reynold stress for the two cases are shown in Figure 13. Generally the agreement between the computed results and the measured data is good. At the region around the vegetation top, the computed results show a sharper variation in the velocity. This is probably due to the fact that the drag characteristics there are affected by the shear at the vegetation top and the local drag coefficient changes. In the vegetation region the velocity is slightly under-predicted for the case with S x = 0.05 m. For the Reynolds stress profiles the peaks occur at the vegetation top. The location of the peak value is accurately predicted (Figure 13(b)). The measured and computed peak values are within 10% difference.
The results clearly demonstrate that for cases of vegetation with the same areal density, the spatial distribution of the vegetation elements can have significant effect on the flow resistance. A reduction in the drag coefficient will reduce the flow resistance and increase the vegetation-layer velocity. The percentage increase in the vegetation-layer velocity is approximately proportional to the square root of the percentage decrease in the drag coefficient. The increase in the vegetation-layer velocity will in turn increase the penetration depth L p and reduce the velocity gradient at the vegetation top, resulting in a moderate increase in the velocity in the clear water region. In the present case with S x = 0.0125 m, the decrease of drag coefficient from 10.0 to 3.1 causes an increase in the depth-average velocity of 43%. For conventional models the drag coefficient is not separately determined and is generally taken to be 2 for blade-type vegetation. For the case with S y = 0.0125 m, the conventional k-model overpredicts the stem-layer velocity significantly and hence underpredicts the hydraulic resistance ( Figure 14).

Conclusions
A DANS model with a k-type turbulence closure has been developed to simulate submerged vegetated flows. The novelty of the model is that the vegetation-generated turbulence dissipation term is related to the areal density of vegetation and the bulk drag coefficient C d . A separate 2D RANS model is used to determine the bulk drag coefficient of vegetation. In contrast to most previous works assuming a constant value of C d for a given areal density of vegetation, the present work shows that C d can have multiple values for a given areal density of vegetation. It slightly increases with increasing the areal density of vegetation (channeling effect) if the transverse spacing is decreased with the longitudinal spacing kept constant. It decreases with increasing areal density of vegetation (sheltering effect) if the longitudinal spacing is increased with the longitudinal spacing kept constant.
Experiments have been conducted to provide additional data for the determination of the turbulence length scale L p used in the DANS model. The turbulence length scale decreases with the drag-vegetation-density parameter f v = C d Nbh v which accounts for the geometrical properties of vegetation. With proper determination of the bulk drag coefficient and turbulence length scale, the DANS model is accurate in predicting submerged vegetated flows.
The results on the variation of drag coefficient with the spatial distribution pattern of vegetation are not exhaustive. The methodology of using a CFD (Computational Fluid Dynamics) approach to compute drag coefficient however is general and has been shown to work properly in the present study. For cases not covered in the present study, such as those in the field, the CFD approach can be used to determine the drag coefficient which can then be input into the DANS model.