Simulation of Water Surface Profiles in Channels with Linear Rigid Emergent Vegetation

1. Introduction

Vegetation along watercourses, on the banks and within floodplain areas, plays a crucial role from both an ecological and environmental perspective. Aquatic plants provide essential nutrients and habitats for animal communities in river ecosystems. Vegetation obstructs sunlight, contributes leaves, and facilitates nutrient exchange. Additionally, it prevents the entry of fertilizers and pollutants into the watercourses. Vegetation also enhances the landscape and holds significant recreational value [1,2,3,4].

Vegetation acts on bank stability since roots increase the soil strength against erosion and reduces the soil water content because of evapotranspiration, so reducing interstitial pressures [5,6,7].

On the other hand, from a hydrodynamic perspective, woody vegetation in the form of trees, shrubs, and bushes occupies part of the river cross-section, increasing roughness and reducing flow velocity. The decreased average velocity reduces riverbed and bank erosion, while simultaneously promoting sediment deposition, leading to a reduction in water cross-sectional area. These changes result in higher water levels and an increased risk of flooding. It is also true that, on the scale of the hydrographic network, the general velocity reduction influences the travel time of water particles, making the peak flow control easier [8,9]. Furthermore, the role of vegetation is increasingly recognized within the broader framework of Nature-Based Solutions (NBSs) for both coastal and marine environments. In these latter contexts, marine vegetation such as seagrass meadows (e.g., Posidonia) acts as a vital vegetative NBS, significantly attenuating wave energy and providing coastal protection [10]. Moreover, vegetation alters the turbulence structure [11,12,13], affects sediment transport [14,15,16,17,18,19], and influences bed morphology [20,21,22]. In conclusion, there is a reciprocal feedback relationship among these various aspects. The resistance provided by vegetation in marshes differs from that in open channels, as flows in marshes typically remain in laminar or transitional states [23].

Historically, in order to maintain water flow, it was common practice for the authorities responsible for watercourse maintenance to remove all or part of the vegetation. However, given the significant role vegetation plays in river ecosystems, the necessity and cost-effectiveness of such interventions have been questioned, particularly due to the high expenses involved. Globally, efforts are now focused on river restoration, re-naturalization, and watershed and watercourse rehabilitation. Effective management requires determining the optimal vegetation density and arrangement to achieve a balance between flood discharge capacity and environmental health. To this end, it is crucial to assess the depths at which a given discharge flows as vegetation density and arrangement vary; in other words, in steady-flow conditions, the water surface profile must be computed. Comparing the effects of different vegetation management strategies allows for the most rational decision-making.

To compute the water surface profile, it is essential to determine the resistance offered by the vegetation. This resistance is a function of various factors, including vegetation flexibility, submergence, foliage, plant morphology, age, mechanical properties, density, and spatial distribution. In the literature, vegetation is typically classified as rigid or flexible, and, depending on water levels, as emergent or submerged.

The objective of this paper is to compute water surface profiles in vegetated channels with rigid, emergent vegetation arranged in a linear pattern, validating the computation results by means of experimental tests. In a recent study, D'Ippolito et al. [24] successfully simulated water surface profiles using a single drag coefficient formula derived for linear arrangements. Building upon this baseline, the present study comprehensively evaluates and compares five different predictive equations. Section 2 briefly outlines the procedure commonly used to compute water surface profiles in channels. In the case of vegetated channels, the need to determine the drag coefficient is emphasized. Additionally, some previously proposed expressions for the drag coefficient are reviewed, particularly for the case of linear vegetation patterns, with a discussion of their applicable ranges. Section 3 presents experimental water surface profiles used to validate the drag coefficient equations considered in Section 2. In Section 4, the experimental profiles are compared with those computed using various drag coefficient equations. This comparison is carried out analytically using statistical parameters and the Taylor diagram. Finally, the physical rationale behind the obtained results is briefly discussed.

2. Methods

It is often necessary to compute the water surface profile in an open channel, in steady-flow conditions for a given discharge, in order to assess the effects of engineering works or modifications to the channel. In this specific study, the reliability of drag coefficient formulas proposed in the literature must be evaluated. The computed water surface profiles, as previously mentioned, are essential for the proper maintenance of watercourses. As is well established, the computation of water surface profiles requires the integration of the following equation using finite-difference methods [25]:

$${\displaystyle \frac{dE}{dx}}=i-J$$

(1)

where x is the coordinate in the flow direction, E is the specific energy, i is the bed slope and J is the slope of the energy grade line. This equation can be integrated using either the direct method or the standard step method. In both cases, the profile computation begins at a control section, which may represent, for example, the critical depth or a depth defined by a hydraulic control structure. The standard step method facilitates the calculation of water depths at specified cross-sections in the case of non-cylindrical channels. In this paper, the method will be employed to compute the water depths at the cross-sections of the laboratory flume where measurements were taken, enabling a direct comparison between the measured and computed depths. In general, J is calculated using the Manning equation. In the literature there are many photos and reference tables for the values of the Manning coefficients in vegetated channels [25,26,27,28,29]. According to some researchers, these indications are more reliable than analytical ones, because direct field observations consider the vegetation heterogeneity. On the other hand, however, photos of vegetated channels are limited and identifying a reference channel is very often difficult or impossible.

For rigid, emergent vegetation with a constant diameter D, the following equation can be applied [30]:

$$J={\displaystyle \frac{C_{D}mD}{\left(1-{\varphi }_{veg}\right)}}{\displaystyle \frac{V^2}{2g}}$$

(2)

where ${\mathrm{C}}_{\mathrm{D}}$ is the drag coefficient, ${\mathsf{\varphi }}_{\mathrm{v}\mathrm{e}\mathrm{g}}$ is the areal concentration of cylinders, $\mathrm{m}$ is the number of vegetation cylinders per unit ground area, V is the mean flow velocity, and g is the gravitational acceleration. Equation (2) underscores the central role of the drag coefficient in the computation of the water surface profile.

In studies involving rigid vegetation, laboratory tests are commonly conducted using wooden or plastic cylinders. This is a physically accurate assumption for representing emergent vegetation (e.g., stiff reeds or tree trunks) where the water flow interacts primarily with the rigid stems and does not submerge the flexible leafage [31]. Furthermore, as recently corroborated by Liu and Wang [32], adopting a simplified layout with purely rigid elements provides an indispensable and fundamental baseline. It allows researchers to properly isolate and investigate the core hydrodynamic mechanisms of canopy drag—such as blockage and sheltering effects—without the confounding variables introduced by the complex bending mechanics of flexible plants. Since the spatial distribution of vegetation significantly influences the drag coefficient, these cylinders are arranged in flumes and studied according to various patterns: linear, staggered, or random (Figure 1).

As it is well known, in the case of a single cylinder, the drag coefficient is a function of the Reynolds number, $Re=VD/\nu$, where $\nu$ is the kinematic viscosity of the fluid [33,34].

When there is a set of cylinders, a mutual influence rises among them, owing to the strong interaction between the wakes as well as between the cylinders and the wakes. The methods used to estimate the drag coefficient can be derived:

• from the force balance per fluid mass e.g., [35];
• based on direct force measurement on a single cylinder e.g., [36], or on a group of cylinders both considering the shear stress e.g., [37], and without considering it e.g., [38];
• on the link with turbulence e.g., [39];
• using computational fluid mechanics e.g., [40];
• employing genetic programming e.g., [41].

For an overview of the various proposed formulas, please refer to D’Ippolito et al. [31], Liu et al. [41] and Gubashi et al. [42].

Below are the formulas related to a linear arrangement of rigid, emergent vegetation that will be used later in the profile simulation.

Wang et al. [43], in a study on the partition of incipient bed shear stress in mobile bed channels, investigated the vegetation drag coefficient for a linear vegetation pattern. The authors conducted 112 tests in a 12-meter-long flume, with a vegetated zone 6 meters in length located in the middle of the flume. The diameter of the cylinders was 6 mm, and the distance between the centers of the cylinders was 5 and 10 cm in the flow direction, and 4, 5, and 6 cm in the transverse direction. The tests were performed without any sediment bed, so the bed shear stress was considered negligible in comparison to the plant shear stress. Based on the experimental data and using a best-fit function approach, the authors proposed the following expression for the drag coefficient:

$$C_D={\displaystyle \frac{90}{{Re}_w^{0.5}}}+4.5{\displaystyle \frac{D}{h}}-0.303ln{\varphi }_{veg}-0.9,$$

(3)

with $h$ the flow depth, ${Re}_w$ a parameter similar to the Reynolds number defined as:

$${Re}_w={\displaystyle \frac{r_v{V}_p}{v}}={\displaystyle \frac{(1-{\varphi }_{veg})V_{p}h}{\nu }},$$

(4)

with $r_v=\left(1-{\varphi }_{veg}\right)h$ representing the original vegetation related hydraulic radius, and $V_p$ the pore velocity, defined as $V_p=V/\left(1-{\varphi }_{veg}\right)$. Equation (3) is considered valid for the range 0.00471 < ${\varphi }_{veg}$ < 0.01413.

D’Ippolito et al. [38] conducted seventy tests with rigid, emergent vegetation, represented by cylinders with diameters of 8 mm and 10 mm. The cylinders were arranged in a square grid, with distances between the centers of the cylinders of 4.24 cm, 8.48 cm, and 12.72 cm. The drag coefficients were calculated from direct drag force measurements, and interpolation of their mean values for each density led to the following logarithmic expression:

$$C_D=0.211ln\left(100{\varphi }_{veg}\right)+0.784,$$

(5)

with R² = 0.92. Equation (5) is considered valid within the following ranges: 0.003 < ${\varphi }_{veg}$ < 0.0436, 0.48% <$i$ < 2.02% and 1000 < $Re$ < 10000.

In a subsequent study D’Ippolito et al. [44] examined twenty-six profiles observed during the seventy tests conducted by D’Ippolito et al. [38] and computed how the drag coefficient varied along each of them. Through multivariate regression analysis, the authors derived the following linear equation

$$C_D=10.83{\varphi }_{veg}+0.2043{\displaystyle \frac{{Re}_D}{1000}}-0.0675{\displaystyle \frac{{Re}_h}{\mathrm{10,000}}}+0.1765,$$

(6)

where ${Re}_h$ is the Reynolds number, ${Re}_h=V_{p}h/\nu$, and ${Re}_D$ is the stem Reynolds number, defined as ${Re}_D=V_{p}D/\nu$, with ${\mathrm{R}}^2$ = 0.83. Equation (6) indicates that $C_D$ increases with ${\varphi }_{veg}$ and ${Re}_D$, while it decreases with ${Re}_h$. This equation was also used to simulate the profiles observed by Wang et al. [30], with satisfactory results when the vegetation density was below 7.3%. Equation (6) is considered valid within the following ranges: 0.0097 < ${\varphi }_{veg}$ < 0.073, 0.0% < i < 2.02% and 2000 < ${Re}_D$ < 7200.

In Equation (6), when computing the free surface profiles, the “effective flow width”, B_e, is used instead of the channel width, B. For a unit length along the streamwise direction, the bed area is B and the area occupied by water is B(1 – ${\varphi }_{veg}$), so the effective flow width is B_e = B(1 – ${\varphi }_{veg}$). As recently demonstrated by D'Ippolito et al. [24], Equation (5) can be successfully extended to higher vegetation density (0.0097 < ${\varphi }_{veg}$ < 0.419) by using the effective flow width.

Based on literature data and the theory by Tanino and Nepf [39], Sonnenwald et al. [45] proposed the following equation:

$$C_D=2\left({\displaystyle \frac{6475D+32}{{Re}_D}}+17D+3.2{\varphi }_{veg}+0.5\right)$$

(7)

where the coefficient of D must have units m^-1 to retain non-dimensionality. The cylinder arrangements representing the vegetation on which Equation (7) is based were square, staggered, or random. Equation (7) is considered valid within the following ranges 0.003 < ${\varphi }_{veg}$ < 0.35 and 12 < ${Re}_D$ < 3838.

Liu et al. [41] employed Genetic Programming (GP) to develop a predictor for the drag coefficient of emergent vegetation directly from published data. This technique automatically identifies all possible relationships between the variables. The authors considered a total of 475 data points from staggered and random arrangements to identify a potential relationship between $C_D$ and three dimensionless parameters - ${{Re}_D,\varphi }_{veg}$ and the blockage factor - without any pre-specified functional forms. The blockage ratio, $\psi =D/L_y$ where $L_y$ is the lateral distance between adjacent cylinders at the same streamwise location, accounts for the difference in drag coefficients between staggered and random distributions at the same density and Reynolds number. The expression for $C_D$ that best balances complexity, accuracy and physical meaning, as proposed by Liu et al. [41] is as follows:

$$C_D={\displaystyle \frac{189}{{Re}_D}}+0.82+{\psi }^2+6.02\lambda .$$

(8)

This equation demonstrated strong predictive capabilities and is valid within the ranges 0.003 < ${\varphi }_{veg}$ < 0.35 and 25 < ${Re}_D$ < 4570.

3. Experimental Data

D’Ippolito et al. [38] performed seventy experimental tests in a variable-slope hydraulic flume, 11.13 m in length, with a PVC bottom, a width of $B=0.382m$, and plexiglass walls 0.21 m in height. The vegetation was simulated using two sets of small wooden circular cylinders with diameters of $D=0.8$ cm and $1.0$cm, placed in the central portions of the flume, with the vegetated section having a length $L$ that varied between approximately 1.8 m and 2.2 m. The cylinders were fixed to two wooden box-structured plates parallel to the flume bottom, positioned above the water surface.

Three different bed slopes were considered: $i=0.48\%,1.35\%$, and $2.02\%$. In all tests, the cylinders were arranged in square meshes with spacing Δx = Δy equal to 4.24 cm (arrangement D9), 8.48 cm (arrangements D4 and D5), and 12.72 cm (arrangement D3). Figure 2 shows arrangements D4, D5 and D9. The experimental tests were conducted with flow rates of approximately Q = 3.0, 5.0, 8.0, 11.0, 14.0 and 16.0 l/s.

The authors obtained 70 values for the drag coefficient through direct measurement of the force F_D on a small group of cylinders, as well as by applying the momentum equation to the control volume containing the aforementioned group of cylinders. These results led to the formulation of Equation (5).

Subsequently, D’Ippolito et al. [44] analyzed 26 water profiles observed during the experimental investigation conducted by D’Ippolito et al. [38]. Many of these profiles exhibited the pattern of an accelerated subcritical flow, type M2, although in some cases, particularly at lower vegetation densities, oscillations were present. In a few instances, an M3 profile, a hydraulic jump, or a gradual transition from supercritical to subcritical flow was observed in the initial part of the reach populated by cylinders or along the length of the flume. Profiles with downstream depths significantly different from the critical depth were also recorded when a sluice gate located at the downstream end of the flume was operated. In these cases, however, the water surface profile corresponded to an accelerated subcritical flow, type M2, often with minimal variations in depth. The water-level profiles were measured using transparent ruler sticks attached to the flume wall, with measurements taken every 4.24 cm. By considering 16 of the 26 type M2 water profiles with linear arrangements D4, D5, and D9, D’Ippolito et al. [44] proposed Equation (6).

In the present study, two additional tests, T30 and T31, were identified as potentially relevant for further analysis, bringing the total number of tests with M2-type water surface profiles to twenty-eight. It is worth noting that 26 of these experimental profiles, alongside the 8 profiles by Wang et al. [30] discussed below, were recently used by D'Ippolito et al. [24] in an investigation to validate solely Equation (5) using the effective flow width. Table 1 summarizes the characteristics of these tests, where each test is identified using the coding system from D’Ippolito et al. [44], with the exception of the newly added tests T30 and T31. The tests labeled T13, T14, and T15 in D’Ippolito et al. [44] were excluded from this analysis, as they correspond to staggered vegetation arrangements. In Table 1, n represents the number of measured depths in the M2 profile section, and h_exp-d denotes the downstream depth. For reference, x = 0 marks the starting point (inlet) of the flow into the flume, and 630 cm corresponds to the abscissa of the initial section of the vegetated zone. For tests T9, T21, and T31, the downstream depths shown in Table 1 were calculated as averages of the last two measurements, due to the oscillations observed in the profiles. In total, Table 1 presents 1052 experimental data points.

Wang et al. [30] conducted eight flume experiments with a flat bed (i = 0) using rigid cylinders with a diameter of D = 8 mm and a length of h_v = 250 mm, representing rigid vegetation with varying densities. The cylinders were arranged in a regular linear configuration similar to that of D’Ippolito et al. [38], but with significantly different vegetation density values. Specifically, the vegetation density ranged from 1% to 41.9%, and the length of the vegetated section varied between 0.52 m and 0.67 m. For all experiments, a steady flow rate of Q = 0.00384 m³ s^-1 was maintained. The flow depths, h(x), for each vegetation density (φ_veg) were recorded using a side-view camera. For reference, x = 0 denotes the starting point (inlet) of the flow into the vegetation zone, where the upstream flow depths, h_exp-u, ranged from 4.7 cm and 21 cm. Some characteristics of the eight tests, including vegetation density (${\varphi }_{veg}$), length of the vegetated section (L), measured upstream water depth (h_exp-u), measured downstream water depth (h_exp-d), critical depth (h_c), and the number of measured depths (n), are provided in Table 2. For tests W7 and W8, the downstream depths presented in Table 2 were calculated as averages of the last four measurements due to the oscillations in the observed profiles. In total, Table 2 includes 444 experimental data points

Wang et al. [30] developed an original procedure to compute the drag coefficient; however, it was not employed in the simulations conducted in this study because its parameters were determined for only a single flow rate.

As observed in the tables above, the vegetation density in the experimental tests conducted by D’Ippolito et al. [38] is similar to that in Wang et al. [30] for only three tests. In the former, the density ranges from 0.78% to 4.36%, while in the latter it ranges from 1.0% to 41.9%. In selecting the equations for simulating the profiles, particular emphasis was placed on vegetation density.

The twenty-eight water surface profiles observed by the authors were simulated using all five Equations (3), (5), (6), (7), and (8) (with some clarifications provided later). The eight water surface profiles observed by Wang et al. [30] were simulated, for vegetation densities equal to or greater than 0.073, using the drag coefficient equations proposed by D’Ippolito et al. [38], Sonnenwald et al. [45], and Liu et al. [41]. The profiles for the experimental tests with vegetation densities less than 0.073 were simulated using all five equations (with some clarifications provided later).

4. Results Analysis

4.1. Results

In the numerical simulations of the profiles, reference is made to the effective flow width, although for the water surface profiles observed by D’Ippolito et al. [44], considering the low vegetation density, its effect is minimal.

The water surface profiles were computed starting from the downstream cross-section. For the M2 profiles presented by the authors, the critical water depth, h_c, was observed at the downstream end of the vegetated section, and is given by:

$$h_c=\sqrt[3]{q^2/\left(g{\left(1-{\varphi }_{veg}\right)}^2\right)},$$

(9)

where $q=Q/B$ is the discharge per unit width, and $g$ is the gravitational acceleration. In some cases, a depth higher than the critical depth was observed (and it was assumed as downstream condition).

Using the standard step method, the twenty-eight experimental water surface profiles presented in Table 1 were simulated. It should be noted that, with reference to arrangements D4 and D5, each test was simulated with all five equations analyzed in this study, as the vegetation density for each experimental test fell within the applicability range of each equation. Tests with arrangement D9 were also analyzed using all five equations. In this case, the vegetation density did not strictly allow the use of the formula by Wang et al. [43], which is valid for vegetation densities between 0.47% and 1.4%; however, its application up to vegetation densities lower than 4.36% provided acceptable results, in line with those obtained using the other formulas. The step size Δx adopted was 1.06 cm.

Figure 3 shows some of the simulated profiles. Specifically, Figure 3a and Figure 3b correspond to arrangement D4, Figure 3c and Figure 3d correspond to arrangement D5, and Figure 3e and Figure 3f, 3g, and 3h correspond to arrangement D9. The other simulated profiles for arrangements D4, D5, and D9 are presented in Sections S1.1, S1.2, and S1.3, respectively, in the Supplementary Information. In Figure 3, as well as in subsequent figures with similar profile comparisons, the scales are intentionally distorted to better highlight the differences between the various simulated profiles. Furthermore, the minimum value on the ordinate axis is slightly lower than the depth measured downstream of the water surface profile.

We first analyze the results obtained for arrangements D4 and D5, followed by those for arrangement D9.

A preliminary analysis of Figure 3a and Figure 3b, 3c, and 3d reveals that for the lowest vegetation densities, specifically those less than φ_veg = 0.0121, the profiles exhibit localized oscillations. These oscillations are, in some cases, concentrated in the central part of the vegetated section, while in others they occur near the terminal end. None of the drag coefficient equations analyzed were able to replicate these oscillations. For all tests conducted with arrangements D4 and D5, a hydraulic jump of varying shape was observed upstream of the M2 profile among the cylinders. Regarding these tests, the plots above display the results of the simulations up to the maximum observed flow depth. However, it is important to note that, since the experimental profiles, as mentioned earlier, exhibit oscillations, and the simulated profiles were computed under the assumption of gradually varied flow, the simulations represent the average trend of the experimental profiles. In cases where the oscillations affected the downstream end of the profile, such as in experimental tests T9, T21, and T31, the average of the last two measured depths was used as the downstream condition, rather than the actual measured depth.

In the computation of the profiles, the drag coefficient equation proposed by D’Ippolito et al. [44] provided overall very good results, while the equation by D’Ippolito et al. [38] yielded acceptable results, with simulated depths slightly lower than the measured ones.

The drag coefficient equations proposed by Wang et al. [43] and Sonnenwald et al. [45] consistently produce simulated depths that are greater than the measured values, as well as those obtained using the other models in this study.

The drag coefficient equation proposed by Liu et al. [41] yields simulated depths that are slightly lower than the experimental ones; however, overall, the simulations are considered good.

We now turn to the results obtained for arrangement D9. A preliminary analysis of the profiles in Figure S1.3 of the Supplementary Information reveals that when the downstream sluice gate is engaged (T5, T11, T20, T23, and T27), the profile exhibits a roughly "linear" trend with a slight bulge in the cross-sections preceding the region where the flow enters the section of the flume without cylinders. In other cases, the profile shows a marked curvature towards the critical depth.

The equations proposed by D’Ippolito et al. [38] and by D’Ippolito et al. [44], when the downstream flow depth is not controlled by the sluice gate, produce simulated profiles that closely match the experimental ones. In the case of "linear" profiles, the formula by D’Ippolito et al. [44] consistently gives lower depths compared to that of D’Ippolito et al. [38], which more accurately simulates the experimental profiles.

The profiles obtained from the Wang et al. [43] formula, which, as previously mentioned, is applied outside the range from which it was derived, provide a generally good estimate of the experimental profiles, particularly for tests T3, T4, T11, T18, T19, T20, T23, T26, T27, T28, and T29. However, it yields profiles with greater depths for the simulations of tests T8, T10, T12, and T17, and with lower depths for tests T1, T2, T5, and T24.

The profiles derived from the formula proposed by Liu et al. [41] appear to most accurately replicate the experimental profiles overall.

The profiles obtained from the drag coefficient formula by Sonnenwald et al. [45] consistently produce higher depths than the experimental profiles, with differences increasing upstream. The results of the simulations using the Sonnenwald et al. ( [45] and Liu et al. [41] formulas are contradictory. Specifically, the profiles based on the Liu’s et al. [41] formula, which apply solely to staggered and random arrangements, are very accurate, while those based on the Sonnenwald et al. [45] formula, which also apply to square arrangements, are less accurate. This discrepancy is likely due to the fact that the data for the square arrangements, from which the Sonnenwald et al. [45] formula was derived, covered only the lowest density and was very limited in number.

Despite these observations, the errors in the various simulations are limited, and overall, all simulations are deemed reasonable.

The profiles identified by Wang et al. [30] were subsequently simulated. The profiles corresponding to vegetation densities of 0.419, 0.291, 0.206, 0.163, and 0.073 were simulated using the drag coefficient equations proposed by D’Ippolito et al. [38], Sonnenwald et al. [45] and Liu et al. [41], even though the vegetation density of 0.419 exceeds the 0.35 density range for which the Sonnenwald et al. [45] and the Liu et al. [41] formulas were developed. The profiles for vegetation densities of 0.041, 0.018, and 0.010 were simulated using all five drag coefficient equations analyzed in this study. The step size, Δx, in the simulation of these profiles was variable and equal to the distance between the cross-sections where the depths were measured. Figure 4 displays both the experimental and simulated profiles, while Section S2.1 of the Supplementary Information presents the water surface profiles for vegetation densities greater than 0.041, which were also simulated using the equations proposed by Wang et al. [43] and by D’Ippolito et al. [44].

The preliminary analysis of the profiles obtained by Wang et al. [30] reveals, despite significant differences in vegetation density, a similar behavior to the profiles obtained by D’Ippolito et al. [44]. Specifically, for vegetation density values above 1.8%, the profiles exhibit a relatively regular shape with a distinct curvature towards the critical depth. For lower densities, the profiles display oscillations. Moreover, for vegetation densities above 16.3%, the results obtained using the formulas proposed by Sonnenwald et al. [45] and Liu et al. [41] are similar.

In the following sections, we first analyze the simulations carried out using the formulas by D’Ippolito et al. [38], Sonnenwald et al. [45], and Liu et al. [41], applied to all experimental tests, and then those conducted using the formulas by D’Ippolito et al. [44] and Wang et al. [43], which were employed in the simulation of experimental tests with vegetation densities less than or equal to 0.041.

The profiles calculated using the drag coefficient equation proposed by D’Ippolito et al. [38] show reasonable agreement with the experimental profiles. The simulations for vegetation densities of 0.419 and 0.073 yield particularly good results. The profile for a vegetation density of 0.291 shows the greatest deviation, with a slight overestimation of the water depths. For lower vegetation densities, the simulations tend to slightly underestimate the water depths.

The formulas proposed by Liu et al. [41] and Sonnenwald et al. [45] produce, for vegetation densities greater than 0.073, simulated depths that are significantly higher than the measured depths, with discrepancies increasing as one moves upstream. As the vegetation density decreases, the differences between the simulated and measured depths gradually diminish. The results are deemed acceptable for vegetation densities of 0.041, 0.018, and 0.010.

Regarding profiles with vegetation densities lower than 4.1%, the drag coefficient formula proposed by D’Ippolito et al. [44] yields slightly smaller depths than the experimental ones, whereas the formula proposed by Wang et al. [43] produces slightly larger depths than those measured.

In summary, for the experimental tests of Wang et al. [30], for vegetation densities greater than or equal to 7.3%, the formula by D’Ippolito et al. [38] provides a good simulation of the experimental profiles. For lower vegetation densities, all formulas offer a reasonable reproduction of the experimental profiles.

4.2. Comparison of the Predicted Water Surface Profiles with the Experimental Ones

As suggested by Pasquino and Gualtieri [46] and by Liu et al. [41], a multiple error criterion should be employed to assess the predictive performance of each drag coefficient (C_D) predictor. For this purpose, we selected the following nine performance indicators:

• Root Mean Squared Error (RMSE);
• Root Mean Squared Error as a percent (RMSE%);
• Mean Absolute Error (MAE);
• Coefficient of Residual Mass (CRM);
• Modelling Efficiency (ME);
• Relative Prediction Error (RPE);
• Mean Bias Deviation (BIAS);
• Correlation Coefficient (R_cc);
• Determination Coefficient (R²).

For further details, refer to Pasquino and Gualtieri [46], Bandyopadhyay et al. [47], Zeng and Huai [48], Vargas-Luna et al. [49]. A complete description of each indicator can be found in the Appendix A, where h_exp,i and h_com,i represent the experimental and computed depths, respectively; ${\overline{h}}_{exp}$ and ${\overline{h}}_{com}$ are their respective averages, and n is the sample size.

For perfect simulation of the observed data, the value of each performance indicator should be as follows: RMSE = 0, RMSE% = 0, MAE = 0, CRM = 0, ME = 1, RPE = 0, BIAS = 0, R_cc = 1, and R² = 1. It is assumed that as the model estimate approaches the observed data, the values of the indicators will approach these ideal values. A positive CRM indicates that the model tends to underestimate the observed values, whereas a negative CRM suggests the model tends to overestimate the observed values.

Furthermore, the model performance was evaluated using the Taylor [50] diagram, which illustrates the standard deviation (STD), centered root mean square difference (RMSD), and correlation coefficient (R_cc) between the measured and estimated profiles.

Figure 5 presents a comparison between the predicted and observed depth values for the 28 water surface profiles observed by the authors, with 0.78% < φ_veg < 4.36%, for each equation used in the computation of C_D.

Starting from the drag coefficient equation proposed by Wang et al. [43], the lower depths of the water surface profiles are overestimated, while the higher depths are underestimated. The profiles obtained using the Sonnenwald et al. [45] equation are systematically overestimated.

The profiles generated from the D’Ippolito et al. [44] equation appear to simulate the experimental profiles less accurately than those derived from the D’Ippolito et al. [38] formula. The drag coefficient equations that provide the best simulations of the experimental profiles are those proposed by D’Ippolito et al. [38] and Liu et al. [41].

Table 3 presents a comparison of the accuracy of the water surface profiles obtained with the various proposed expressions for C_D.

For each indicator, the best value is highlighted in bold. As it is evident, the equation that most accurately simulates the profiles is that of Liu et al. [41], with values closely approximating the desired ones and an RMSE of 3.54%. Following this is the equation of D’Ippolito et al. [44], which yields comparable indices,

Regarding the water surface profiles derived from the equation proposed by Sonnenwald et al. [45], it is noteworthy that the CRM value is negative (with an absolute value higher than any other), indicating, as observed in the previous graph analysis, that the model tends to overestimate the experimental values, with an RMSE of 9.43%.

The Wang et al. [43] equation results in depths with an RMSE of 7.4%.

The negative CRM value of -0.03 can be attributed, as observed in the graph in Figure 5a, to the overestimation of smaller depths, particularly those less than 10 cm, which occurs more consistently than the underestimation of larger depths, approximately 10 cm. This outcome is not surprising, however, since, as previously noted, the Wang et al. [43] equation was applied to vegetation density values outside the range for which it was originally calibrated.

The statistical parameters indicate that the D’Ippolito et al. [44] formula simulates the experimental profiles more accurately than the D’Ippolito et al. [38] formula, which contrasts with the findings from the graph in Figure 5. Specifically, the RMSE for the former is 3.92, while for the latter, it is 5.11.

The high values of the correlation coefficients and the coefficients of determination should also be emphasized for all simulations, as they indicate a strong relationship between the experimental and calculated values, primarily due to their consistent trend. However, these two coefficients do not provide a proper assessment of the validity of the different formulas. They will, nonetheless, be utilized in Taylor's [50] diagram in subsequent analyses.

Figure 6 presents a comparison between the predicted and observed depth values corresponding to the water surface profiles from Wang et al. [30], for each equation used in the computation of the drag coefficient (C_D). Section S2.2 of the Supplementary Information displays the same figure, but with the drag coefficient equations from both Wang et al. [43] and D’Ippolito et al. [44] applied to all vegetation densities.

The results obtained using the D’Ippolito et al. [38] equation appear to be quite satisfactory, although a closer inspection of the figure reveals a slight overestimation in some of the tests.

As shown in Figure 6, the equations of Sonnenwald et al. [45] and Liu et al. [41] tend to significantly overestimate the experimental depths.

The graphs for the D’Ippolito et al. [44] and Wang et al. [43] equations display fewer data points, as only experimental tests W6, W7, and W8 were simulated. However, they indicate that the profiles are reasonably simulated, with a slight underestimation for the D’Ippolito et al. [44] formula and a slight overestimation for the Wang et al. [43] formula.

A detailed analysis of Figure 6b, Figure 6d and Figure 6e, which are based on a significantly smaller set of experimental data compared to the simulations for the D’Ippolito et al. [44] experimental test, reveals that the points are aligned along different segments. This is due to the fact that, as one moves upstream, the depths either increase or decrease and progressively deviate from the measured values, exhibiting a variation that is approximately linear.

Table 4 presents a comparison of the accuracy of the water surface profiles obtained using the different proposed expressions for the drag coefficient (C_D).

The results regarding the accuracy of the simulations conducted using the different formulas, as shown in Table 4, cannot be directly compared because they correspond to different amounts of data, specifically 172 and 444 points. Comparisons can, however, be made within each group separately. The first group includes the formulas of D’Ippolito et al. [38], Sonnenwald et al. [45], and Liu et al. [41], while the second group consists of the formulas of D’Ippolito et al. [44] and Wang et al. [43].

Regarding the first group, it is evident, as previously noted in the analysis of Fig. 4, that the expressions of Sonnenwald et al. [45] and Liu et al. [41] significantly overestimate the depths, resulting in RMSE values of 35.81% and 34.15%, respectively. This is further confirmed by the CRM values of -0.22 and -0.19, respectively.

Overall, the formula of D’Ippolito et al. [38] provides the best simulation of the experimental profiles, with an RMSE of 5.58%.

A comparison of the results obtained using the drag coefficient equations proposed by D’Ippolito et al. [44] and Wang et al. [43] for the test simulation with a vegetation density of less than 0.041 shows that the former yields slightly better results, with an RMSE coefficient of 5.29%, compared to 6.16% for the latter.

Taylor diagrams (Figure 7) are commonly used to assess the performance of multiple predictors, as they succinctly summarize the degree of agreement between observed and simulated data. The diagram is presented in a polar format, where the azimuth angle represents the correlation coefficient (R_cc) between the observed and simulated values. The radial distances from the origin to the points representing the observed and simulated values are proportional to their respective standard deviations. The point representing the standard deviation of the observed values is plotted along the abscissa. In Fig. 7, the dashed lines represent the distance from the reference point and indicate the root mean square difference (RMSD), which approaches zero as the observed and simulated models become more similar.

Figure 7a presents the Taylor diagram for the simulations of the water surface profiles by D’Ippolito et al. [44]. It illustrates that the formula proposed by Liu et al. [41] provides the best simulation of the experimental profiles. However, the formulas proposed by D’Ippolito et al. [38] and D’Ippolito et al. [44] seem to yield comparable results.

The Taylor diagram for the Wang et al. [30] tests, shown in Figure 7b, was developed using only the simulation results obtained with the drag coefficients proposed by D’Ippolito et al. [38], Sonnenwald et al. [45], and Liu et al. [41]. From this diagram, it is evident that the drag coefficient formula by D’Ippolito et al. [38] delivers the best results. Overall, the equations of Liu et al. [41] and Sonnenwald et al. [45] provide comparable results, as also observed in the analysis of the profiles and statistical parameters. In the case of vegetation densities of less than 4%, the equation by Liu et al. [41] yielded the best simulation of the experimental profiles.

Section S2.2 of the Supplementary Information presents Figure 7b, but with the drag coefficient equations from Wang et al. [43] and D’Ippolito et al. [44] applied to all vegetation densities.

In summary, it can be concluded that the equation proposed by D’Ippolito et al. [38] provides overall good results across the entire range of analyzed vegetation densities. This equation has the advantage of being independent of the Reynolds number. It is recommended for use in the vegetation density range ${0.01\le \varphi }_{veg}\le 0.41$ and for Re_D > 1000, and it enables the computation of type M2 water surface profiles.

However, in the authors' opinion, it would be prudent to verify these results under different bed slopes and flow rates. To the best of the authors’ knowledge, there are no other available data in the literature on tests involving rigid, emergent vegetation arranged in a linear distribution, apart from those of Zhang et al. [51]. These profiles, however, pertain to regularly arranged natural reeds, for which the density is reported in terms of reed trunks. This, however, is not particularly significant, as each trunk often possesses several branches.

5. Discussion

The analysis of the simulated profiles and statistical parameters revealed that the effectiveness of the drag coefficient formulas is highly dependent on the vegetation density range.

The formula proposed by Wang et al. [43] provides good estimates for densities lower than 0.041, but it tends to overestimate the minor depths and underestimate the major depths. This outcome is expected, given that the equation was applied outside the range of vegetation density for which it was originally calibrated. From a hydrodynamic perspective, this inaccuracy at higher densities is linked to the blockage effect. The equation was calibrated for a very narrow and low density range (${0.0047\le \varphi }_{veg}\le 0.014$), where individual cylinder wakes behave almost independently. At higher densities, the physical blockage constricts the effective flow width, accelerating the pore velocity and drastically altering the pressure distribution around the stems. These non-linear hydrodynamic phenomena cannot be captured by a formula calibrated only on sparse vegetation.

Contradictory results were obtained when evaluating the formulas by Sonnenwald et al. [45] and Liu et al. [41]. Both significantly overestimate the experimental depths for vegetation densities greater than 0.073. As demonstrated by Etminan et al. [52], in linear arrangements - such as those investigated in this study - the cylinders are aligned parallel to the flow direction, creating a strong sheltering effect where upstream cylinders reduce the flow velocity and modify the turbulent wake directly impacting the downstream cylinders. Conversely, formulas derived from staggered or random arrangements (such as Sonnenwald et al. [45] and Liu et al. [41]) fail because the flow is forced into a more tortuous path. Because downstream cylinders are positioned directly in the accelerated flow paths between the upstream elements, the longitudinal sheltering is effectively interrupted and the actual incident velocity increases [32]. This geometric shift makes the blockage effect dominant, severely constricting the cross-sectional area and naturally leading to a much greater flow resistance [52]. However, for very low densities (less than 4%), the distance between elements is large enough that the difference in flow resistance due to the spatial arrangement is reduced, allowing the equation by Liu et al. [40] to provide optimal simulations.

Overall, the formula by D’Ippolito et al. [38] demonstrated the highest robustness, providing very good results over the entire range of density analysed (${0.0078\le \varphi }_{veg}\le 0.419$), with the advantage of being independent of the Reynolds number. This highlights a crucial physical aspect: in the considered subcritical accelerated flows with rigid emergent cylinders, the flow regime is fully turbulent (Re > 1000). In such conditions, the inertial form drag strongly dominates over the viscous friction drag. Therefore, the total drag coefficient becomes predominantly a function of the geometric blockage and vegetation density, rather than viscous scale effects, fully explaining the robustness of this Reynolds-independent equation.

6. Conclusions

The aim of this paper was to assess the ability of several equations proposed in the literature for the drag coefficient in order to simulate water surface profiles in vegetated channels with a linear distribution of vegetation. Indeed, most of the aforementioned equations were derived under uniform flow conditions, and there are only a few experimental studies in the literature that focus on profile simulation, despite such profiles being common in natural environments.

Three of the formulas used in this study (Wang et al. [43]; D’Ippolito et al. [38]; D’Ippolito et al. [44]) were derived from experimental tests with a linear distribution of vegetation, one (Sonnenwald et al. [45]) from data based on linear, staggered, and random distributions combined, and the final one (Liu et al. [41]) from staggered and random arrangement data.

The experimental profiles were obtained from the studies by Wang et al. [30] and D’Ippolito et al. [44], along with two additional profiles observed by the authors. The water surface profile simulations were conducted using the standard step method. The experimentally measured profiles had vegetation densities that did not allow all of the aforementioned formulas to be applied for the simulation of each profile. For vegetation densities below 4.36%, the simulations indicated that the formulas derived from linear distributions generally provided good results, except for the formula proposed by Wang et al. [43], which was applied outside the range for which it was originally derived. The formulas proposed by Liu et al. [41] and Sonnenwald et al. [45] yielded good results for vegetation densities less than 4.4%, but overestimated the depths for higher densities. This is due to the fact that the resistance of staggered vegetation is significantly higher than that of linear vegetation. Overall, the profiles obtained using the formula proposed by D’Ippolito et al. [38] provided very good results.

The results presented above were further validated through the analysis of statistical parameters and the Taylor [50] diagram. This was made possible by the use of the standard step method for the computation of water surface profiles, which provided the calculated depths precisely at the abscissae where the measurements were taken. Finally, the underlying physical explanations for the results obtained were analyzed.

In conclusion, the results obtained using the formula proposed by D’Ippolito et al. [38] across the entire range of vegetation density variability ($0.0078\le {\varphi }_{veg}\le 0.419$) are generally satisfactory. The proposed equation is independent of the Reynolds number, and thus unaffected by water viscosity, making it easily applicable to similar arrangements within the specified range of densities and for ${Re}_D\ge 1000$.

However, it is believed that the aforementioned formula should be further validated through new experimental tests conducted at higher vegetation densities and with bed slopes other than zero.