Friction factor for turbulent open channel flow covered by vegetation

The need for operational models describing the friction factor f in streams remains undisputed given its utility across a plethora of hydrological and hydraulic applications concerned with shallow inertial flows. For small-scale roughness elements uniformly covering the wetted parameter of a wide channel, the Darcy-Weisbach f = 8(u*/Ub)2 is widely used at very high Reynolds numbers, where u* is friction velocity related to the surface kinematic stress, Ub = Q/A is bulk velocity, Q is flow rate, and A is cross-sectional area orthogonal to the flow direction. In natural streams, the presence of vegetation introduces additional complications to quantifying f, the subject of the present work. Turbulent flow through vegetation are characterized by a number of coherent vortical structures: (i) von Karman vortex streets in the lower layers of vegetated canopies, (ii) Kelvin-Helmholtz as well as attached eddies near the vegetation top, and (iii) attached eddies well above the vegetated layer. These vortical structures govern the canonical mixing lengths for momentum transfer and their influence on f is to be derived. The main novelty is that the friction factor of vegetated flow can be expressed as fv = 4Cd(Uv/Ub)2 where Uv is the spatially averaged velocity within the canopy volume, and Cd is a local drag coefficient per unit frontal area derived to include the aforemontioned layer-wise effects of vortical structures within and above the canopy along with key vegetation properties. The proposed expression is compared with a number of empirical relations derived for vegetation under emergent and submerged conditions as well as numerous data sets covering a wide range of canopy morphology, densities, and rigidity. It is envisaged that the proposed formulation be imminently employed in eco-hydraulics where the interaction between flow and vegetation is being sought.

Since its inception by Darcy (in 1857) and Weisbach (in 1845), the now called Darcy-Weisbach equation for determining frictional losses in open channels (and pipes) is considered 'standard' provided its associated friction factor coefficient = gR S h f is the friction velocity related to the kinematic 'bed stress' , g is the gravitational acceleration, R h is the hydraulic radius, S f is the frictional slope or energy grade-line that converges to the bed-slope S o for uniform flow, is the bulk or time and area-averaged velocity determined from the flow rate Q and cross-sectional area = A Bh w (assuming rectangular section) orthogonal to the flow direction, B is the channel width and h w is the water depth approximating R h when  B h w . The definition of f can be combined with the estimate of u * to yield = . S fU gR 8 (1) In classical hydraulics, the work of Moody, Nikuradse and many others established f to vary with two dimensionless quantities: the relative roughness r/R h and a bulk Reynolds , where ν is the kinematic viscosity of water. This expression is generally accepted in pipe-and open channel-flows above small-scale roughness elements where  r h / 1 by the so-called Strickler scaling ∼ f r R ( / ) h 1/31-3 . An immediate consequence of this result is that when r is a priori known, f and subsequently S f can be determined from Eq. 1.
Operationally, such an expression for S f can be used to mathematically close the combined continuity and unsteady shallow water flow equations (i.e., the Saint-Venant) to predict h w and U b in a plethora of hydrological and hydraulics applications. Example applications include overland flow from bare to vegetated patches, dam-break and similar shallow inertial flows, flash-flood runoff in ephemeral streams, tsunamis on coastal plains, to name a few [4][5][6][7][8][9][10][11] . It is now recognized that such naive view of f cannot be juxtaposition to streams covered by large roughness values ( > . r h / 01 w ) such as shallow flow over gravel beds [12][13][14][15][16][17] or vegetation [18][19][20] and frames the scope of the work here. In fact, a number of authors are calling for the abandonment of such an approach altogether in such situations 21,22 .
The main theoretical novelty is to arrive at a new expression for f whose generic form resembles the conventional = * f u U 8( / ) b 2 for steady-uniform flow through or above dense canopies. Specifically, it is shown that under certain simplifying assumptions, the canopy-related f (hereafter referred to as f v ) is given by where U v is the spatially averaged velocity within the canopy volume, and C d is a local drag coefficient per unit frontal area derived to include the layer-wise effects of vortical structures within and above the canopy along with key vegetation properties. For this expression to be used in practice, an estimate of U v /U b and C d are required. The derivation of Eq. 2 and estimates of C d and U v /U b are first presented followed by a comparison between model predictions and measurements of f v . Comparisons with other widely used models of f v are also featured demonstrating that the proposed formulation appears superior to all prior formulations across a wide-range of experiments and canopy-flow configurations.
Background. To progress on the description of f v for canopy flows, a number of studies have been conducted to explore connections between the shape of the mean velocity profile u(z) and its depth-integrated value defined as b w h 0 w in wide rectangular channels, where z is the distance from the channel bottom [23][24][25][26][27][28][29][30] . These studies naturally bridge the dominant vortical structures 31,32 to u(z) 12,33-35 and subsequently to U b and f v . Other studies have focused on the force balance and stresses acting on the vegetation elements as well as R h to arrive at u * and f v [36][37][38] . Various resistance models for flow within and above vegetated elements were proposed that used bulk flow measurements or combination of models and measurements to determine f 9,39-56 . The emerging picture from all these studies is that key vortical structures impact the shape of u(z) across various canopy layers and in the vegetation free layer for submerged vegetation. The relation between these vortical structures and u(z) is now reviewed.
Review of the mean velocity profile within and above canopies. For stationary and planar-homogeneous flow in a wide channel covered by a densely vegetated canopy, the flow region h w can be decomposed into three zones according to the dominant sizes of the vortices in each zone as shown in Fig. 1.
Zone I forms near the channel bottom where von Karman vortex streets dominate the energetics of turbulence in the vertical direction. The mean velocity in zone I can be approximated by a constant given as 41,57,58 = u gS C mD where C d is a local drag coefficient as before, m is the number of rods or vegetation stems per unit ground area, and D is the frontal width of vegetation elements. www.nature.com/scientificreports www.nature.com/scientificreports/ Zone II spans the vegetation top and is dominated by attached eddies to a zero-plane displacement as well as mixing-layer eddies (though those types of eddies do not co-exist in space but both impact time-averaged statistics). When explained by the mixing layer analogy, the mean velocity is given by where u v,top is the mean velocity at the top of the vegetation elements, h v is the vegetation height, L mix is a characteristic energetic eddy size generated from Kelvin-Helmholtz instabilities 34 .
Zone III resembles a canonical turbulent boundary layer above a rough ground surface and is commonly represented by for the surface layer (as before) whose thickness is h s , z o is the momentum roughness length (linked to r or h v ), d is the zero-plane displacement height and is related to the penetration depth, κ is the von Karman constant, and w v (submerged vegetation case). Due to the presence of vegetation, the addition of zones I and II introduce major distortions to the expected u(z) shape when compared to a typical rough boundary (i.e. zone III) where f v primarily varies with r/h w and Re b,h as before. When  r h / 1 w (in zone III), r may be linked to z 0 using the scaling relation ∼ r z o 1/6 discussed elsewhere 12 . The aforementioned distortions to u(z) and the possibility of deriving a general expression for f v to be used in operational models for emergent and submerged vegetation is the main goal of the work here.

Derivation of a Friction Factor for Vegetated Flow
The goal of the derivation here is to arrive at a new formula for f v whose generic form resembles the conventional = ⁎ f u U 8( / ) b 2 and is given by This definition necessitates estimates of the friction velocity for the vegetated flow = where S v denotes the energy grade-line slope caused by the presence of vegetation, R v denotes a new vegetation-related hydraulic radius, and U b is, as before, the bulk flow velocity averaged over the entire cross-sectional area but adjusted for the finite porosity of the vegetation medium. Section 2.1 derives the vegetation-related energy slope S v and Section 2.2 derives the vegetation-related hydraulic radius R v . The outcome is then compared to a wide range of published data sets and prior formulations derived from fitting to other published data sets.
Derivation of vegetation-related energy slope. Throughout, it is assumed that the flow is steady and uniform occurring within a wide rectangular open channel characterized by width B and bed slope S o . The case considered here is when vegetation is sufficiently dense so that the overall friction factor is mainly due to f v not bed and side-wall friction and the frictional slope is S v . The vegetation is further assumed to be cylindrical with height h v , diameter D, and density m defined by the number of cylinders per unit area as before. The vegetation can be in one of two-states: submerged (α < 1) or emergent (α > 1) as characterized by the degree of submergence α = h h / v w . For a control volume above a unit bed area extending from the channel bed to the water surface, the mean momentum balance in the streamwise direction reduces to a force-balance between (i) the gravitational contribution of the water weight along the streamwise direction for a unit ground area (=F w ), and (ii) the drag force per unit ground area (=F v ) caused by the presence of vegetation stems resisting the flow. The F w is given by where ρ is the water density, V w is the volume of flowing water per unit ground area determined from the entire domain volume reduced by the volume occupied by the vegetation elements. The αφ for emergent vegetation (i.e. α > 1), where φ is the volume fraction of the vegetation. When ignoring wall stresses relative to the drag force imposed by the dense vegetation and assuming a quadratic-law for F v yields where A v is frontal area of vegetation given by is a 'local spatially-averaged' non-constant drag coefficient that must account for the effects of the vegetation on the flow as discussed elsewhere 55 [59][60][61] . (b) The pore velocity in VL calculated based on a spatially-averaged value www.nature.com/scientificreports www.nature.com/scientificreports/ where the effective width used for the pore velocity is The constricted velocity in a classical staggered array 55 , where L s,stag is the spacing distance defined in Etminan et al. ' work 55

. (d) The separation velocity
, where k p is a kinetic energy of a subgrid scale (SGS) for a Smagorinsky model 55 . Due to different vegetation arrangements used in the experiments analyzed here (Section 3 and 4), the spatially-averaged pore velocity was adopted (i.e. . Detailed analysis shows that this choice of U ref along with a vegetation-related hydraulic radius resulting in a new Reynolds number perform reasonably when compared to many experiments (discussed later). Hence, from the force balance = F F w v , the vegetation-related friction slope is given by where C d and U v,p remain to be determined.

Derivation of a vegetation-related hydraulic radius R v .
Many studies treat flow resistance of a vegetated canopy as equivalent to a 'bed stress' assuming energy losses occur due to a 'vegetated bed' . This equivalence is convenient as it leads to a hydraulic radius related to flow depth h w or surface layer depth h s 9,41,44,66 without considering the details of the vegetation configuration such as frontal area associated with φ, m and D. However, the flow resistance by vegetation is mainly dominated by form drag and is directly linked to the frontal area of the rods. The hydraulic radius may be revised to account for these vegetation-related features. Here, a vegetation-related hydraulic radius R v is proposed by extending the work of Cheng 36 from emergent condition to submerged cases, where the vegetation-related hydraulic radius R v is now defined by the ratio of the whole flow volume V w per unit ground area to the vegetation-fluid contact A resistance per unit ground area and can be expressed as From Eq. (9), the vegetation resistance is caused by its form drag, which is a function of the frontal vegetated area A v per unit ground area. Hence, the vegetation-related hydraulic radius is Since both V w and A v are defined per unit ground area, the ground area cancels in the definition of R v .
Derivation of the friction formula for flow through submerged and emergent vegetation. The friction factor caused by the presence of vegetation is now obtained by substituting Eqs (10 and 12) into Eq. (7) to yield This expression is similar to a prior result given by 9 to be the velocity difference between the vegetation zone U v and the free water zone U s . Equation (14) can be arranged to read where C d is a local drag coefficient per unit frontal area (as shown next). More relevant here is that Eq. (13) allows for multiple effects to be conveniently included in a dimensionless C d (instead of a dimensional L c ) and U v /U b . Before doing so, it is to be noted that when submergence and the conventional expression for friction factor is recovered from Eq. (13). To show explicitly similarities and differences between Eqs (13 and 14), the force balance for submerged conditions is rewritten per unit ground area as where C d is a local drag coefficient per unit frontal area as before. Equation (14) further assumes the hydraulic radius to be the flow depth, or However, when www.nature.com/scientificreports www.nature.com/scientificreports/ using a vegetation-related hydraulic radius αφ is recovered where C d remains a local drag coefficient per unit frontal area (as before). Naturally, when setting (13) is recovered. It can be surmised that different R h definitions yield expressions for f that appear to be different. However, not withstanding these apparent differences, the C d remains a local drag coefficient per unit frontal area in all of them.
Because C d is to be related to a local Reynolds number, several possibilities are now reviewed for defining a local Reynolds number in the presence of vegetation. Etminan et al. 55 compared Reynolds numbers using various characteristic velocity scales (U v,b , U v,p , U v,c and U v,s ) with a fixed characteristic length scale D giving Re v,b , Re v,p , Re v,c and Re v,s . The aforementioned work then showed that typical C d formulation for a single cylinder can be employed to predict a drag coefficient for staggered vegetation when using U v,c as the reference velocity to form and resulting in a best-fit expression applicable for < < Re 0 6000 v c , given as where suffix 'EGL' denotes the first letter of the surname of each author in Etminan et al. 55 . Here, vegetation effects are introduced in both velocity and length scales by selecting the pore velocity (U v or U b ) and the vegetation-related hydraulic radius R v for defining a Reynolds number. For emergent canopies, only the canopy layer exists and the Reynolds number is labeled as Re v,v (in the suffix, the first letter 'v' is vegetation layer, and second letter 'v' denotes the characteristic length in the Reynolds number to be the vegetation-related hydraulic radius) given as where R v is the vegetation-related hydraulic radius as before, given as For a submerged vegetation, the Reynolds number Re b,v can also be defined as (in the suffix the first letter 'b' is bulk flow, and second letter 'v' denotes the characteristic length in the Reynolds number to be vegetation-related hydraulic radius) where R v is, as before, the vegetation-related hydraulic radius now given as

Friction Factor for Flow Through Emergent Vegetation
Cheng's result 36 for emergent vegetated flow is recovered. In the aforementioned study, f was shown to be a monotonically decreasing function of Reynolds number Re v,v consistent with several studies 29,39,59,62,63,67,68 . The relation between f and Reynolds number is analyzed using several experiments described next. www.nature.com/scientificreports www.nature.com/scientificreports/ A local C d,local -Re function in steady non-uniform flow was derived and shown to be parabolic in shape (i.e. not monotonic) within the vegetation zone for a given vegetation density. Here, the focus is on the averaged drag for the entire vegetation zone, which can be expressed as where C d,local is a local drag coefficient that varies in the streamwise direction along the vegetation zone, the normalized distance is = + x x L / veg along the streamwise direction, and L veg is the streamwise length of the vegetation zone. A summary of the data used here is given in Table 1 and the details are shown in Supplementary Information (Table S1). Also, further details about the experimental setup can be found elsewhere 51 .
Two other published data sets from Ishikawa et al. 59 and Tanino et al. 62 are used here and are briefly summarized in Table 1.

Experiments conducted by Ishikawa et al.
The first data set was generated from a 15 m long, 0.3 m wide, and 0.3 m deep channel using two types of steel cylinders with diameters 0.004 m and 0.0064 m, and a height of 0.2 m as described elsewhere 59 . The cylinder spacing was uniform in each case (=0.0632 m and 0.0316 m, respectively). The flow was nearly uniform and the steel canopy array was arranged in staggered configuration. The hydraulic parameters required here are given in Supplementary Information (Table S1). Tanino Table S1) were tracked from their best-fitting curve to yield an approximate expression where suffix 'tn' denotes the expression of Tanino et al. 62 , T 1 and T 2 are obtained from linear regression. The Re v,d in their formulation was based on a stem Reynolds number given as

Experiments conducted by
The linkage between vegetation-array and stem related Reynolds number is Here, a simplified expression that summarizes all the aforementioned expressions may be expressed as a function of Re v,v (desirable for the purposes of the work here) as = .
+ . . C Re 0 819 58 5 This summary expression is shown in Fig. 2, which an be re-arranged to be related to the more commonly reported Re v,d as

Determination of depth-averaged velocity in the vegetation and surface layers.
In the vegetation layer (VL), the depth-averaged velocity can be derived based on the force balance between vegetation drag and flow gravity in streamwise direction For surface layer (SL), the depth-averaged velocity can be determined by the linkage between the vegetation-layer velocity U v and the bulk velocity U b . The bulk velocity is defined by the total discharge Q to the effective cross-sectional area and given as (pore velocity for bulk flow) From the continuity equation,  www.nature.com/scientificreports www.nature.com/scientificreports/ Arranging the above equation, the bulk velocity is given by α αφ Rigid vegetation  www.nature.com/scientificreports www.nature.com/scientificreports/ then the depth-averaged velocity of surface layer gives αφ α Methods to determine friction factor in the submerged cases. The measured friction f v,measure can be obtained from U v (using Eq. 34), U b (using Eq. 35) and C d (using Eqs 30 or 31) when setting = f 2 . All data points and their variation with the Reynolds number Re b,v or Re b,d are given in Fig. 4. The results suggest no obvious trends between measured f v and estimated Re b,v or Re b,d .
The absence of a unique relation between f v and Reynolds number requires further inquiry. The key to calculating f v is U v /U b . Previous studies 9 showed that f v is linked to two dimensionless groups: α = h h / v w and η = h L / v c . As mentioned before, the adjustment canopy length scale is defined as The aforementioned study 9 did show a nonlinear increase in f v with increasing α at a given η, and a nonlinear increase in f v with increasing η at a given α. Based on these results, we propose a 3-parameter mathematical function to describe this behavior without focusing on the detailed mean velocity profile in each zone. This 3-parameter function is given as  www.nature.com/scientificreports www.nature.com/scientificreports/ where p 1 , p 2 , p 3 are constants to be determined from regression analysis. The corresponding friction f v can now be expressed as Friction factor derived from previous studies. Some studies report expressions for the velocity within the vegetated zone, which can also be used to obtain a friction factor. A summary of these expressions are briefly reviewed. Stone et al. 74 estimated the bulk velocity as Baptist et al. 41 report a bulk velocity as where C d is the bed-related Chezy coefficient (=60 m 1/2 s −1 for smooth bed).
Huthoff et al. 28 proposed a Yang et al. 82 showed a bulk velocity given as   www.nature.com/scientificreports www.nature.com/scientificreports/ Katul et al. 44 proposed an analytical model for U b linked to three canonical length scales (L c , h v , h w ) and is given as where U KCL , U KSL are layer-averaged velocity for the canopy and surface layers, respectively, given as . min m D (0 135 , 0 33) is a momentum absorption coefficient, and z 0 is given by With these mean velocity profile formulations, friction factor can be estimated. The bulk velocity of the vegetation layer is determined using Eq. (34) and the velocity across h w is calculated based on Eqs (43 to 48).
Comparison between measured and modeled friction factor using published bulk velocities. The different methods are now used to predict f v and are then compared with measured f v in Fig. 6. To assess the agreement between prediction and measurement, R cc is computed and reported in Table 3. Figure 6 suggests that the present model performs 'no worse' than other (more elaborate) methods. In many cases, it even outperforms prior methods.

Discussion and Conclusion
2 encodes much of the recent developments about vegetation effects on local drag (e.g. blockage and sheltering, weak Reynolds number dependence, etc) and their up-scaled contribution on Many studies found that C d in a vegetated array differs from C d,iso . When only focusing on variations in vegetation density, C d appears to increase 62,68 and then decrease 60,84 with increasing vegetation density 55 . In the present study, we go beyond vegetation density and introduce all the key vegetation effects in a Reynolds number formed by U v and R v . These velocity and length scales contain the vegetation density φ and appear to collapse the experiments onto a single curve.
To accommodate the mechanics of sheltering, delayed separation and blockage effects arising from an array of cylinders instead of isolated cylinders, C d of Eq. (31) can be used. Sheltering effect indicates that some elements were located in the wake region of the upstream elements 85 , resulting in a lower velocity than their upstream counterparts and generate a lower drag compared with the isolated cylinder case. Delayed separation can be explained by the enhancement of the mean separation angle that is larger than that for the isolated cylinder, resulting in a decreasing drag coefficient (compared with the isolated cylinder). Both sheltering and delayed separation tend to reduce drag when compared to the isolated cylinder case. For the blockage effect, which leads to a local increase in the drag coefficient, it can by explained by two main factors 55 , one is that the velocity between cylinders is enhanced by the presence of the vegetation in flow. The other factor is reduced wake pressure 86 .
The C d,iso is a local drag coefficient acting on a single stem per unit frontal area without considering the interaction among elements in a vegetation array. To illustrate the role of sheltering, delayed separation and blockage on C d , a 'bifurcation' type Reynolds number + Re v d , may be introduced. When < + Re Re v d v d , , , the viscous boundary layer formed around cylinders creates a slow-moving flow zone with a path smaller than the spacing between adjacent stem. This effect leads to an equivalent local drag coefficient by the array of cylinders to be larger than the one associated with isolated cylinders (i.e. a blockage effect). However, when > + Re Re www.nature.com/scientificreports www.nature.com/scientificreports/   Table 3. Comparison between measured and modeled friction factor by different formulations.
www.nature.com/scientificreports www.nature.com/scientificreports/ Karman streets that grow in size and 'fill' the space between the stems. This effect is akin to a decreasing drag coefficient when compared with C d,iso (i.e. sheltering effect). Blockage and sheltering effects can be distinguished by calculating When > E 1, blockage dominates and when < E 1, sheltering dominates C d 51 . Using Eqs (31 and 52) with different vegetation concentration (φ = 0.01, 0.05, 0.1, 0.3 and 0.5), blockage and sheltering can be delineated for Re v,d ranging from 10 1 to 10 5 in Fig. 7. As expected, C d increases with increasing φ for small φ. The threshold Reynolds number + Re v d , also increases with increasing φ. For example, at ≈ + Re 4000 v d , , separation between blockage and sheltering occurs for φ = .
0 01, when , , the C d for isolated stems are comparable but not identical to their array counterpart.
To conclude, the expression of friction factor proposed here does accommodate blockage or sheltering through a local drag coefficient C d (Eq. 31) and any distortions to the shape of the mean velocity profile (Eq. 41). Evaluated using numerous data sets covering a wide range of canopy morphology, densities and rigidity, the friction factor for emergent vegetation is shown to be proportional to the drag coefficient. This finding shows why the friction factor monotonically decreases with increasing Reynolds number for emergent vegetation. For submerged vegetation, the friction factor is shown to vary with submergence and adjustment length scale. Also, this variation is not monotonic with increasing Reynolds number.