A genetic programming-based model for drag coefficient of emergent vegetation in open channel flows

Highlights

• A new Ergun-like model for estimating drag coefficient of emergent vegetation is developed.

• Machine learning technology is introduced into the study of drag coefficient.

• The dependence of the two parameters in Ergun equation on vegetation characteristics are determined by genetic programming.

Abstract

The estimation of drag exerted by vegetation is of great interest because of its importance in assessing the impact of vegetation on the hydrodynamic processes in aquatic environments. In the current research, genetic programming (GP), a machine learning (ML) technique based on natural selection, was adopted to search for a robust relationship between the bulk drag coefficient (C_d) for arrays of rigid circular cylinders representing emergent vegetation with blockage ratio (ψ), vegetation density (λ) and pore Reynolds number (Re_p) based on published data. We utilize a data set covering a wide range of each parameter involved to cover all possible dependencies. A new predictor, which shares the same form with the Ergun-derived formula, was obtained without any pre-specified forms before searching. The dependence of the two parameters in Ergun equation on vegetation characteristics was also estimated by GP. This new C_d predictor for emergent vegetation with a relatively concise form exhibits a considerable improvement in terms of prediction ability relative to existing predictors.

Introduction

Aquatic vegetation plays an important role in river ecosystems and greatly influences hydrodynamics processes within aquatic environments because of the additional resistance exerted by emergent and submerged vegetation (Etminan et al., 2017; Green, 2006; Huai et al., 2012; Nepf, 1999; Stoesser et al., 2010). The additional drag caused by vegetation may considerably reduce the flow velocity of the river, raise the water level and increase the risk of flooding. Moreover, by modifying the velocity profiles and turbulence structures in vegetated channels, the aquatic vegetation helps to reduce sediment erosion and suspended transport. Thus, a comprehensive understanding of vegetation resistance is essential to the assessment of the impact of vegetation on these hydrodynamic processes in aquatic environments. Many researchers have a strong interest in the estimation of the drag generated by vegetation because of its great significance in ecological and environmental engineering (Mattis et al., 2012; Tang et al., 2014).

In the past few decades, research on vegetation drag has mostly been carried out in laboratories. Given the complex geometry of real vegetation in the field, aquatic plants are usually modeled as an array of rigid circular cylinders with uniform diameter, d. This simplified geometry of real vegetation has been widely recognized by most researchers (Hamidifar et al., 2016; Li et al., 2018; Li and Zeng, 2009;Liu et al., 2008). The effect of canopies on the flow is mainly determined by the additional resistance, which can be quantified by the drag coefficient, defined as:

$$C_d=\frac{F_D}{\frac{1}{2}\rho A{U}_{ref}^2}$$

where F_Dis the vegetation-induced drag, ρ is the water density, A is the frontal area of the vegetation stems (In the case of an isolated circular cylinder, A=dL, L is the length of the cylinder immersed in water.), and U_ref is a reference flow velocity. Usually the mean velocity approaching the cylinders is used as the reference velocity, which is equal to the upstream free flow velocity in the case of an isolated cylinder. For emergent vegetation, given the impact of vegetation distribution, using upstream free flow velocity (herein referred as “bulk flow velocity”, U) as the reference velocity may be reasonable only for very sparse vegetation. However, the mean pore velocity U_p, $U_p=U/\left(1-\lambda \right)=Q/\left(Wh\left(1-\lambda \right)\right)$, is closer to the flow velocity that is actually approaching the vegetation stems than the bulk flow velocity, because the time-averaged longitudinal flow velocity is almost uniform in the vertical direction within emergent vegetation (Liu et al., 2008). Here, W is the width of the channel; Q is the channel flow rate; λ is the dimensionless vegetation density ($\lambda =\pi n{d}^2/4$, where d is cylinder diameter and n is the number of cylinders per unit bed area). Therefore, in most studies, the mean pore velocity U_p is typically used as the reference flow velocity (Kothyari et al., 2009; Tanino, 2012; Tanino and Nepf, 2008).

Usually, the drag coefficient C_d for emergent vegetation is obtained from laboratory experiments either by directly measuring the drag forces using load cells (strain gauges) or by equating driving forces with resistance caused by cylinders. The first method includes two forms. In one form, the drag force generated by the hole cylinder array is directly measured with a drag plate (Tinoco and Cowen, 2013; see Fig. 1 (b)), which can better reflect the bulk drag characteristics of the array. Obviously, the resistance measured by this method includes vegetation and bed resistance. In the numerical simulation results of Stoesser et al. (2010), bed shear stress only accounts for approximately 0.82%–6% of the total flow resistance; thus, bed resistance can be completely ignored. However, the disadvantage of this method is the relatively high requirement for experimental equipment. A simpler alternative method is to measure the drag force generated by a single cylinder in the array (Ishikawa et al., 2000; Kothyari et al., 2009; van Rooijen et al., 2018; see Fig. 1 (a)). That is, we only need to select one cylinder in the area where the flow is fully developed to measure its drag and use it to represent the drag characteristic of the entire cylinder array. More details about these two different methods can be found in Tinoco and Cowen (2013) and Kothyari et al. (2009), respectively. When the drag force is measured, C_d can be obtained from Eq. (1).

When the experimental conditions do not allow the direct measurement of the drag force, C_d can be derived from the force balance per fluid mass within the vegetated channel, which can be expressed as (Nepf, 1999):

$$\left(1-\lambda \right)C_b{U}_p^2+\frac{1}{2}{C}_{d}ah{U}_p^2=gh\left(1-\lambda \right)\frac{\partial \eta }{\partial x}$$

where the first term on the left is the resistance exerted by the channel bed, with a bed friction coefficient, C_b, the second term on the left is the drag force contributed by the vegetation stems; λ is vegetation density, a is frontal facing area (For circular cylinders, $\lambda =\pi ad/4$ and a = nd.), h is flow depth, g is acceleration due to gravity, U_p is mean pore velocity, x is the streamwise direction, and η is water surface elevation. As mentioned above, the bed shear stress only accounts for 0.82%–6% of the total flow resistance, therefore, the first term on the left of Eq. (2) can be omitted. Then, C_d can be written as:

$$C_d=\frac{2gS\left(1-\lambda \right)}{a{U}_p^2}$$

where S is the energy slope, $S=\partial \eta /\partial x$. Obviously, when the uniform flow condition is reached in the experiment, S is just equal to the bottom slope of the channel (e.g., Cheng and Nguyen, 2011; James et al., 2001, 2004; Kim and Stoesser, 2011). For non-uniform flow, when the channel bed is horizontal, S is equal to ∂h/∂x (e.g., Ferreira et al., 2009; Meftah and Mossa, 2013; Nepf, 1999; Tanino and Nepf, 2008; Tanino, 2012).

Moreover, drag coefficient can also be obtained directly from high resolution numerical simulations, such as wall-resolving Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS) (Chang et al., 2018; Etminan et al., 2017; Nicolle and Eames, 2011; Stoesser et al., 2010).

However, physical measurements and numerical simulation approaches are often unavailable or infeasible for the practical application of ecological and environmental engineering. Therefore, C_d must be estimated. Researchers have approached the study of estimating drag coefficient for emergent vegetation from different angles: pure theoretical derivation (Meftah and Mossa, 2013), regression analysis on the basis of a pre-given form obtained from theoretical derivation (Sonnenwald et al., 2017,2018; Tanino and Nepf, 2008), and multi-parameter regression analysis without any physical basis (Kothyari et al., 2009).

The drag coefficient predictor for an isolated cylinder is well established which can be expressed as (White, 1991):

$$C_d=1+10R{e}^{-2/3}$$

where the cylinder Reynolds number is $Re=\rho Ud/\mu$, μ is dynamic viscosity. Note that for an isolated cylinder, λ=0, i.e., U=U_p. However, the flow structure in a multi-cylinder array is much more complicated than the flow structure around a single isolated cylinder because of the strong wake-wake interaction and wake-cylinder interaction, especially when the vegetation density is not very small. Thus, Eq. (4) may not be directly used to estimate the drag coefficient of the multi-cylinder array. Based on Eq. (4), different C_d predictors for cylinder array can be obtained by defining Reynolds number and drag coefficient with various velocity and length scales (i.e., reference velocity and reference length). In emergent canopies, as mentioned above, previous researches have chosen the reference velocity scale to be the mean pore velocity, U_p. However, through LES, Etminan et al. (2017) found that the “constriction cross-section velocity”, U_c, is the velocity scale that actually governs wake pressure and thus vegetation resistance. Therefore, Etminan et al. (2017) and van Rooijen et al. (2018) expressed the canopy drag coefficient, using the new velocity scale as:

$$C_{d,c}=1+10R{e}_c^{-2/3}$$

where C_{d, c} and Re_c are the canopy drag coefficient and cylinder Reynolds number based on U_c. $U_c=\left(1-\lambda \right)U_p/\left(1-\sqrt{4\lambda /\left(\beta \pi \right)}\right)=U/\left(1-d/L_y\right)$, β is the ratio of the lateral distance between adjacent stems at the same streamwise location, L_y, to the distance between two rows of vegetation stems in the streamwise direction, L_x (as is shown in Fig. 5).

Further, Cheng (2013) re-fitted the empirical formula of the C_d for an isolated circular cylinder using the available experimental data. Subsequently, a pseudofluid model is developed to define a non-dimensional cylinder diameter, $d_{*}={\left(2C_{d}R{e_p}^2/\pi \right)}^{1/3}$. With this new dimensionless cylinder diameter, Cheng (2013) extended the drag coefficient formula that was used for an isolated cylinder to the cylinder array. Finally, the drag coefficient formula for cylinder array can be written as:

$$C_d=\left(\begin{array}{c}11{\left(\frac{\left(1+S\right)R{e}_p}{1+80\lambda }\right)}^{-0.75}+0.9\left[1-\exp \left(-\frac{1000\left(1+80\lambda \right)}{\left(1+S\right)R{e}_p}\right)\right]\hfill \\ +1.2\left[1-\exp \left(-{\left(\frac{\left(1+S\right)R{e}_p}{4500\left(1+80\lambda \right)}\right)}^{0.7}\right)\right]\hfill \end{array}\right)\frac{1+S}{1-\lambda }$$

where Re_p is cylinder Reynolds number based on the mean pore velocity U_p and cylinder diameter d, i.e., $R{e}_p=\rho U_{p}d/\mu$.

As mentioned above, by using different velocity and length scales to define Reynolds number and drag coefficient, different C_d predictors can be obtained. Etminan et al. (2017) and van Rooijen et al. (2018) redefined Reynolds number and C_d using a new velocity scale, U_c, and obtained Eq. (5). Similarly, Cheng and Nguyen (2011) proposed a new length scale, vegetation-related hydraulic radius, r_v, to characterize vegetated channels. The new Reynolds number defined with the new length scale is written as:

$$R{e}_v=\frac{\rho U_p{r}_v}{\mu }=\frac{\pi \left(1-\lambda \right)}{4\lambda }R{e}_p$$

where $r_v=\pi \left(1-\lambda \right)d/\left(4\lambda \right)$.With this new Reynolds number, Re_v, Cheng and Nguyen (2011) proposed a best-fit function to empirically describe the relationship of C_d and Re_v based on existing experimental data of C_d. The best-fit function can be expressed as:

$$C_d=\frac{50}{R{e_v}^{0.43}}+0.7\left[1-\exp \left(-\frac{R{e}_v}{15000}\right)\right]$$

Ergun (1952) proposed an expression for pressure drop in packed columns, and Tanino and Nepf (2008) related this expression to drag coefficient giving:

$$C_d=2\left(\frac{{\alpha }_0}{R{e}_p}+{\alpha }_1\right)$$

where α₀ and α₁ are coefficients describing the viscous contribution that arises from the viscous shear stress on the cylinder surface and the inertial contribution that arises from the pressure loss in the cylinder wake, respectively. Sonnenwald et al. (2018) re-parametrized α₀ and α₁ in terms of λ and d based on least-squares curve fitting using available experimental data. Then, they obtained a new function for estimating C_d:

$$C_d=2\left(\frac{6475d+32}{R{e}_p}+17d+3.2\lambda +0.5\right)$$

where the coefficients of the d terms must have units m⁻¹ to retain dimensionless. Tanino and Nepf (2008), Tinoco and Cowen (2013), and Sonnenwald et al. (2017) obtained different expressions with similar forms, which are not listed here.

Kothyari et al. (2009), using their experimental data, proposed a C_d predictor by multi-parameter regression analysis without considering any physical background written as:

$$C_d=1.53\left[1+0.45\ln \left(1+100\lambda \right)\right]R{e_p}^{-3/50}$$

Obviously, all drag coefficient predictors mentioned above were developed with data-driven methods, such as regression analysis or curve-fitting with pre-specified forms. However, they were obtained with a relatively limited variable space and not physically sound, so they can only perform well in some certain cases (as discussed in Section 3.2). When extrapolated to conditions outside the sampled parameter space, it often produces absurd drag coefficient values or makes the predictors meaningless (e.g., Eqs. (8) and (10)). Generally, researchers analyze the dependence of drag force on flow parameters and vegetation geometry and then subjectively suggest the variables that may affect C_d, and even further, subjectively suggest an analytical expression, i.e., giving a pre-specified form of C_d predictor (Sonnenwald e al., 2017, 2018). This means that the predictive ability of the model is closely related to the reasonableness of the subjective “suggestion”. For instance, Sonnenwald et al. (2018) pre-specified that the C_d predictor has the Ergun-like form and considered that α₀ was linearly related to d, and α₁ was linearly related to d and λ. Then a drag coefficient model with mixed dimensional and non-dimensional parameters was obtained, however, it may lead to ridiculous C_d values (considering a large cylinder with 1 m diameter), which is dangerous. Moreover, several predictors developed by traditional multi-parameter regression methods are very complicated and not physically interpretable (e.g., Eqs. (6), ((8), and (10)), which is not conducive to our further understanding of the influence mechanism of flow variables and vegetation geometry on drag force.

However, a more advanced data-driven method, machine learning (ML) technique, enables the exploration of all possible relationships among variables in complicated phenomena on the basis of a provided data set. In present research, we utilize genetic programming (GP), an ML technique based on natural selection to develop a predictor of drag coefficient for emergent vegetation directly from the experimental data that have been published before. Compared with the conventional data-driven methods mentioned above, GP does not require researchers to make relatively subjective analysis of the data set. It is based on the principle of natural selection to automatically find all possible relations between the variables and only retains the solutions with higher accuracy in the evolving process, instead of needing researchers to subjectively suggest a possible expression of drag coefficient on flow variables and vegetation features. At the same time, the GP algorithm can take all variables into the development of a predictor, while leaving the task of identifying related variables to the program (Wang et al., 2017). In contrast, as mentioned above, the conventional regression methods usually rely on the researchers’ subjective “suggesting”. Nevertheless, the current understanding of the influence mechanism of drag force is not yet sufficient. As a result, these subjective suggestions will lead to a poor predictor. While the GP algorithm can well avoid the subjective suggesting from researchers. Furthermore, since the time used to analysis the data set manually is greatly reduced in GP, more time can be devoted to explore the physical background of C_d (Tinoco et al., 2015; Wang et al., 2017). Thus, the solution given by GP is simple in form and physically sound, thus having good performance.

In this study, a well-documented and widely tested GP software, Eureqa (Schmidt and Lipson, 2013), which has been successfully applied to the studies of vegetated flow (Huai et al., 2018; Tinoco et al., 2015; Wang et al., 2017), is used to find a new predictor without any pre-specified forms (in contrast to previous data-driven methods which pre-specify a form and proceed to find the optimal relationship through regression analysis. e.g., Sonnenwald et al., 2017, 2018; Tanino and Nepf, 2008; Tinoco and Cowen, 2013).