*Article* **Deciphering Morphological Changes in a Sinuous River System by Higher-Order Velocity Moments**

## **Jyotismita Taye 1, Jyotirmoy Barman 1, Bimlesh Kumar 1,\* and Giuseppe Oliveto <sup>2</sup>**


Received: 29 January 2020; Accepted: 9 March 2020; Published: 11 March 2020

**Abstract:** Bank erosion in a sinuous alluvial channel is a continuous phenomenon resulting in bank instability and migration of sediment. In this study, flume experiments were conducted in a sinuous channel to investigate its morphological changes and hydrodynamics. High-order velocity fluctuation moments are analyzed at outer and inner banks to explain the morphological variation in a sinuous river channel. The variance of streamwise velocity fluctuations on both banks of the sinuous channel follows a logarithmic law from a particular depth. In the outer bend region, the magnitude of velocity fluctuation moment is significant, indicating erosion. The trend of velocity fluctuation at higher even-order moments is similar to the variance of streamwise velocity fluctuations where the outer bend magnitude is greater than the inner bend. The premultiplied probability density functions (PDFs) and the flatness factor show greater magnitude in the outer bend of the channel as compared to the inner bend.

**Keywords:** sinuous channel; velocity fluctuations; river bend erosion; structure function

## **1. Introduction**

The study of turbulence in a sinuous channel is a complex topic as compared to a straight channel. There is uniformity in a straight channel in both banks; however, in a sinuous channel, there are continuous erosion and deposition processes at outer and inner banks, respectively. Turbulence studies in a sinuous channel have been discussed quite elaborately over the past 2–3 decades. Researchers such as Rozovski˘ı [1], Anwar [2], de Vriend and Geldof [3], Blanckaert and Graf [4], Booij [5], Sukhodolov and Kaschtschejewa [6], and Engel and Rhoads [7] have studied the turbulent characteristics of both infield and laboratory and forwarded theories and articles regarding the uniqueness of flow behaviour. Several researchers (e.g., [8,9]) have also explained the relations between turbulence characteristics and erosion at the outer bank. Blanckaert [10] stated that the Reynolds shear stress (RSS) in a sinuous bend indicates the presence of helical flow, which contributes to the erosion and deposition processes. Esfahani and Keshavarzi [11] studied the bursting process by octant analysis. Their study was based on two models of 17◦ and 30◦ bend angle. They found that the effect of river bends on flow characteristics and bursting events is inversely proportional to the curvature of the bend. The velocity distribution at channel bends is conflicting. In sinuous bends, the maximum velocity is not found always towards the outer bend because of the velocity redistribution [3,12]. de Vriend and Geldof [3] mentioned that there is a shift of maximum velocity in the inner bend before the flow enters the bend apex. This paper aims to address that the mean velocity may not be always higher towards the outer bend. Thus, the paper highlights the erosion and deposition patterns in a sinuous channel by analyzing the high-order streamwise velocity fluctuations.

Previous studies have testified the morphological changes in a sinuous channel [13–15]. In sinuous channel, the shifting of centerline and migration of sediment is a continuous process. Experiments conducted with erodible beds in a meandering channel have reported bed formations as point bars, mid-channel bars, free bars, and sand bars [16–19]. Xu and Bai [12] in their experiments with erodible bed and fixed walls, observed depositional bars at the inner bank and pools at the outer bank. Zhang et al. [20] observed the flood events in a step-pool channel. Scour appeared in the pool at threshold condition, and after which the scour depth increased linearly until the step collapsed. Riverbanks are vulnerable to erosion, and recent studies have been done in riverbanks of Langat River in Malaysia [21] and Parlung Tsangpo River in China [22]. Under natural conditions in rivers, the flow turbulence governs the bed deformation and channel migration [23]. Small-scale physical models of sinuous channels were prepared successfully in the laboratory [24–27], although it is very challenging to simulate a dynamic realistic sinuous channel. However, to find a solution to a specific problem to understand the fundamentals of river processes, experimental channels are quite useful to explore the mechanisms associated. In this study, scale effects were neglected, and a laboratory study was conducted to assess flow turbulence effects.

Many researchers have worked towards high-order velocity moments in the turbulent boundary layer to acquire the information on flow behaviour. Meneveau and Marusic [28] studied the higher-order moments of velocity fluctuations at different Reynolds numbers and presented a generalized logarithmic law. de Silva et al. [29] studied the statistical properties of wall turbulence using higher-order moments of streamwise velocity fluctuations. They found that at higher Reynolds numbers, the even-order moments follow the logarithmic behaviour. Sharma and Kumar [30] studied the higher-order structure functions for seepage-affected channels. They found that higher-order moments of velocity fluctuations increase with the application of downward seepage. Another parameter called the flatness factor or Kurtosis helps in observing the distribution of the fluctuating velocity. Meneveau and Marusic [28] in their study estimated the Reynolds number dependence on sub-Gaussian statistics. The flatness factor evaluated in the inertial region is below the Gaussian value 3 and found no noticeable dependence of the Reynolds number on sub-Gaussianity.

The previous studies present the significance of the high-order velocity moments and the correlation of the turbulent flow with sediment transport. Previously, outer bend erosion was explained by using various turbulent parameters such as Reynolds shear stress, turbulent intensities, and turbulent kinetic energy. However, the previous studies had not addressed the scope to understand the fluvial morphology in association with the higher-order velocity fluctuations. Here, we have considered only the streamwise velocity to analyze the higher-order velocity fluctuations. Though the velocity field is three-dimensional in our channel and has an impact on the morphology, the study has been limited to acknowledge specifically the scour and deposition at the outer and inner bend, respectively. In this paper, firstly, we have analyzed the mean streamwise velocity and the morphology of the sinuous channel. The morphological behaviour around the bend apex (i.e., scour and deposition) needed more insights, for which the mean velocity distribution is not sufficient to make a statement on the bed deformation. Concerning this, we have applied the high-order structure function to understand the scour and deposition in a sinuous bend.

#### **2. Experimental Methods and Program**

The experiments were conducted in a glass-sided recirculating flume of length 17.2 m and width 1 m (Figure 1). In the main channel, we constructed a rectangular sinuous channel with rigid sides of length 5.64 m and width 0.3 m. The centerline of the sinuous channel follows the sine-generated function forwarded by Leopold and Langbein [31]. The function is expressed as

$$\theta = \theta\_0 \cos \left( 2\pi \frac{\text{m}}{\text{L}} \right) \tag{1}$$

where θ<sup>o</sup> is the maximum angle the channel makes with the downvalley axis, θ is the angle at distance m measured along the centerline of the channel, and L is the arc length (length along the centerline of the channel between two repeating points). The channel centerline at the crossover made an angle of 65◦ with the horizontal, known as the deflection angle. The channel sinuosity was calculated as σ = L/λ [32], where L and λ are the arc length and wavelength of the meander, respectively (Figure 1b). The wavelength was calculated as λ ≈ 2πB [33], where B = 0.3 m is the width of the sinuous channel. The calculated sinuosity and wavelength were equal to 1.25 and 1.88 m, respectively. At the entry and exit of the sinuous channel, guide vanes were provided to direct the flow into and out from the sinuous channel smoothly. The flow was administered into the main channel from an overhead tank with the help of a control valve. The water first falls into an inlet tank (2.8 m length, 1.5 m width, and 1.5 m deep) at the upstream of the flume and it gradually enters the main channel. A tailgate at the downstream of the flume manages the flow depth during experiments. The main channel discharge was measured using a rectangular notch located at the downstream end of the channel. The experiments were conducted on a uniform river sand bed of median diameter d50 = 1.1 mm. The flow discharge attained was equal to 0.0156 m3/s with a flow depth of 0.117 m. The flow depth was measured along the centerline of the channel. In our study, the channel achieved the Reynolds number and Froude number as 57, 731 and 0.40, respectively. Therefore, we have maintained a subcritical flow condition (Froude number < 1) and turbulent flow condition (Reynolds number > 10, 000). We have physically modeled a sinuous channel in a laboratory flume. Scale effects may arise in physical modeling of sediment transport processes when all the forces in the model and the real field river flows are nonidentical. As gravity is the primary driving force in open channel flows, we tried to achieve a Froude number similar to field conditions. In the field study by Engel and Rhoads [23], the onsite Froude number is about 0.3 (subcritical flow), and the Reynolds number ranges from 292, 698 to 397, 254 (turbulent flow). The Froude number observed in our tests is very close to the field conditions. Therefore, scale effects due to this may be negligible. However, the Reynolds number of the laboratory flow is greater than 50, 000. Most of the river flows are turbulent and in the hydraulic rough regime, where losses are independent of the Reynolds number. Therefore, the Reynolds number of the laboratory flows was greater than 50, 000 and the shear Reynolds number (R<sup>∗</sup> = 75) was achieved greater than 70 so that the laboratory flow was in the fully turbulent hydraulic rough regime to better account for the losses [24].

The instantaneous flow velocities u, v, and w in three directions X (streamwise), Y (transverse), and Z (vertical) were measured by a velocimeter. The instantaneous velocities are equal to the sum of mean velocities (u, v, and w) and fluctuating velocities ( u , v , and w ) in the form u = u + u , v = v + v , and w = w + w . The experimental investigation of recording the velocity and morphological changes were performed along the second bend of the sinuous channel. We used an acoustic doppler velocimeter (ADV) by Nortek®(vectrino+ 4-beam down-looking probe, Nortek AS, 1351 Rud, Norway) to record the velocity data at five locations (1, 2, 3, 4, and 5) of the bend at section r − r, s − s, and t − t (Figure 1c). The velocimeter works on the principle of Doppler effect and has four downward-looking probes with the sampling volume located 0.05 m beneath the central transmitter. Readings were taken at the five locations throughout the flow depth. At each point of measurement, the velocimeter recorded 12, 000 samples for 120 s (sampling rate 100 Hz). For higher-order statistical analysis, Schwarz et al. [34] recommended that at least 10, 000 samples should be collected.

The uncertainty associated with the ADV measurements was tested by taking 17 pulses for a duration of 120 s recorded at near-bed depth z ∼ 3 mm, where z is upward and is positive in the vertical direction (Table 1). The data collected from the ADV contains spikes, and therefore it should be filtered. The data were filtered using the acceleration threshold method [35]. During the measurement, the signal-to-noise ratio (SNR) and correlation were greater than 15 dB and 60%–70%, respectively. The correlations were reduced by ± 5% approximately near the channel bed [36]. The spikes were filtered such that the despiked data satisfies the Kolmogorov's 5/3 law in the inertial subrange (Figure 2), where the acceleration threshold value ranged from 1 to 1.5 by trial and error [37].

**Figure 1.** Schematic diagram of the experimental setup showing: (**a**) Side view of the experimental flume; (**b**) plan view of the experimental setup; (**c**) section r − r, s − s, and t − t where velocities were measured at locations 1, 2, 3, 4, and 5; (**d**) sections from "a" to "q" where ultrasonic ranging system (URS) readings were made to track morphological changes.

**Table 1.** Uncertainty test of acoustic doppler velocimeter (ADV) data.\*

**Figure 2.** Power spectra Fuu(f) cm2/s of unfiltered and filtered velocity time-series at outer and inner bends. Fuu(f) is the velocity power spectra of the streamwise velocity u, which is a function of frequency f (in Hz).

To examine the morphological changes along the bend, we used the ultrasonic ranging system (URS). It is a SeaTek®(1.0 cm diameter model, Seatek SPa diesels, Lombardy, Italy,) manufactured instrument consisting of eight transducers. In this instrument, the transducer acts as both transmitter and receiver. The transducer first transmits a pulse of 10-microsecond duration, and then this pulse travels through the water and reflects off a target. The reflected signal travels back to the transducer and is detected back by the electronics. The system has an accuracy of ±0.2 mm. The uncertainty associated with the URS measurements was evaluated by taking 16 sample recordings at the center of the bend (Table 2). The URS measured the bed elevation along the bend in 17 sections (Figure 1d). The transducers were mounted on a trolley and were aligned horizontally in a single line to track the changes. The URS tracks the distance (in centimeters) between the channel bed and the water surface.

**Table 2.** Uncertainty test for ultrasonic ranging system (URS) measurements.


## **3. Results and Discussions**

The velocity profile in the outer (location 1) and inner (location 5) bend of the section s − s of the sinuous channel is shown in Figure 3. The mean velocities u, v, and w in streamwise, transverse, and vertical directions are calculated as:

$$\overline{\mathbf{u}} = \frac{1}{\mathbf{n}} \sum\_{\mathbf{i}=1}^{\mathbf{n}} \mathbf{u}\_{\mathbf{i}} \tag{2}$$

$$\overline{\mathbf{v}} = \frac{1}{\mathbf{n}} \sum\_{\mathbf{i}=1}^{\mathbf{n}} \mathbf{v}\_{\mathbf{i}} \tag{3}$$

$$\overline{\mathbf{w}} = \frac{1}{\mathbf{n}} \sum\_{i=1}^{n} \mathbf{w}\_i \tag{4}$$

where n is the total sample number. In Figure 3, we observe that the magnitude of inner velocity is more when compared to the outer bend of the channel throughout the flow depth. This result is similar to that found by the authors of the papers [3,12,38] who claimed that the velocity is greater in the inner bend than outer because of velocity redistribution. Rozovski˘ı [1] considered the logarithmic distribution for the streamwise velocity profile expressed as:

$$\frac{\mathbf{u}}{\overline{\mathbf{U}}} = 1 + \frac{\sqrt{\mathbf{g}}}{\kappa \mathbf{C}} (1 + \ln \overline{\mathbf{z}}) \tag{5}$$

where U is the depth-averaged streamwise velocity, κ is the von Kármán constant, C is the Chézy coefficient, and z = z/h (z is the height to the point of measurement above channel bed, and h is the flow depth). The experimental profiles of the streamwise velocity in our study show an acceptable correlation with Equation (5).

**Figure 3.** Streamwise velocity profile at the outer (location 1) and inner bends (location 2) of the sinuous channel.

Contour plots show the mean velocity u (m/s) distribution across three sections r − r, s − s, and t − t (Figure 4). The mean velocity is maximal towards the inner bend. In section t − t, the velocity is distributed throughout the channel width. This finding reveals that the mean velocity is not always greater towards the outer bend. Due to the inward skewing, the main velocity may take longer time to reach the outer bend, and therefore the maximum velocity is redistributed mostly towards the inner bend. Shams et al. [38] observed higher streamwise velocity towards the inner boundary in their study on a physical and laboratory-scale model.

**Figure 4.** Contour plots of streamwise velocity at the bend cross-sections (**a**) r − r,(**b**) s − s, and (**c**) t − t of the sinuous channel.

The morphological changes are also analyzed along a bend of the sinuous channel. The variations in an alluvial bed are visible in the channel. The contour plots of the morphological changes were represented using the Surfer®(Golden Software, Colorado, US) [39] at different time intervals (Figure 5). After the desired discharge was achieved in the channel, the morphology readings were taken at intervals of 2, 6, and 10 h Readings were taken up to 12 h, and after which no significant changes were noticed in the bend. It can be seen that the outer bend experiences erosion, which is increasing over time [40,41]. With the acquired flow discharge, visible transport of the sediment took place. The flow interaction in bends allows the sediment to move in a transverse or radial direction (perpendicular to the direction of flow) across the bend. This motion of sediment is due to the established secondary currents in bends. Previous studies [8,10,15] have already focused on the average turbulent parameters such as the bed shear stresses, Reynolds stresses, and secondary currents to explain the scour and deposition in bends. Here, we have focused on the high-order turbulence characteristics and how they affect the morphological processes in a sinuous river. Investigation at 2 h shows the initial development

of variation in the bed along the outer and inner bend. With time (6 and 10 h), the scour depth at the outer bend is prominent.

**Figure 5.** Morphological changes along the second bend of the sinuous channel after 2, 6, and 10 h. This bend was selected because it was unaffected by the entry and exit conditions.

The cross-sectional bed elevation across the bend apex (section i) is shown in Figure 6. Considering the initial level as the datum, the outer bend experiences erosion, which has increased with time. At the outer bend of section i, the depth of scour after 2 h run was found to be 4.35 cm. After 10 h run, the scour depth was estimated to be 7.99 cm. The extreme lower point in the vertical axis is the maximum depth the scour has reached. There is sediment below this point, and thus the scour has not touched the rigid bed of the channel. The flow characteristic and the channel planform play a significant role in the development of morphological changes in a channel. However, the observed morphological changes in the sinuous channel do not conform to the velocity distribution given in Figures 3 and 4. Hence, this conflicting behaviour might be explained by higher-order velocity moments.

**Figure 6.** Cross-sectional morphological changes across the bend apex (section i) after 2, 4, 6, 8, and 10 h.

The velocity profile in both the outer and inner banks varies, as the pressure forces on the banks are different. The logarithmic law of mean velocity profile in the inertial region (i.e., inner flow zone) is given as:

$$\frac{\langle \mathbf{u} \rangle}{\mathbf{u}\_{\star}} = \left\langle \mathbf{u}^{+} \right\rangle = \kappa^{-1} \ln \left( \frac{\mathbf{z} \mathbf{u}\_{\star}}{\mathbf{y}} \right) + \mathbf{B} \tag{6}$$

where u<sup>∗</sup> [=(τ◦/ρ) 0.5] is the shear velocity and τ<sup>0</sup> (bed shear stress) is found out by using the TKE (turbulent kinetic energy) method, υ is the kinematic viscosity of water, z is the wall distance (i.e., distance from the bed), κ is the von Kármán constant, and B is the constant. The value of u<sup>∗</sup> and τ<sup>0</sup> for the channel are 0.034 m/s and 1.2 N/m2, respectively.

Many previous efforts by various researchers were put forward to understand the erosion and deposition behaviour in the outer and inner bend of a sinuous channel. In this section, this behaviour of a sinuous channel will be looked upon from the perspective of structure function. Here, we establish a relation between higher-order velocity moments with the erosion and deposition across sinuous bend. The erosion is a result of the transport of sediment particles. From the granular perspective, the motion of a sediment particle depends upon the balance of drag force exerted by the fluid flow and submerged weight of the particle. The classical turbulence parameter associated with transport is the average bed shear stress τo. When τ<sup>o</sup> exceeds τ<sup>c</sup> (critical bed shear stress), erosion is expected. In

our experiments, the sand bed was in motion. Moreover, the tractive stress depends upon the fluid velocity above the particle <sup>τ</sup><sup>o</sup> <sup>∼</sup> u2. Therefore, the behaviour of instantaneous component of velocity fluctuations becomes significant to erosion, as the instantaneous drag force is directly dependent on it. Hence, we have analyzed higher moments of turbulence in the inner and outer bend.

Studies by various researchers [42–44] have found a logarithmic nature of velocity fluctuations, which is given as:

$$
\left\langle \left( \mathbf{u'}^{+} \right)^{2} \right\rangle = \mathbf{B}\_{1} - \mathbf{A}\_{1} \ln(\mathbf{z}/\delta) \tag{7}
$$

where <sup>u</sup><sup>+</sup> = (u−<sup>u</sup> ) <sup>u</sup><sup>∗</sup> is the nondimensional fluctuating component of streamwise velocity and <sup>δ</sup> is the boundary layer thickness. The velocity fluctuating moments raised to the pth root follows logarithmic nature as proposed by Meneveau and Marusic [28]:

$$\left\langle \left( \mathbf{u}^{\prime +} \right)^{2\mathbf{p}} \right\rangle^{1/\mathbf{p}} = \mathbf{B}\_{\mathbf{p}} - \mathbf{A}\_{\mathbf{p}} \ln(\mathbf{z}/\delta) = \mathbf{D}\_{\mathbf{p}}(\mathbf{R}\mathbf{e}\_{\star}) - \mathbf{A}\_{\mathbf{p}} \ln \mathbf{z}^{+} \tag{8}$$

where z<sup>+</sup> = zu<sup>∗</sup> ν is the distance to the wall and Dp <sup>=</sup> Bp <sup>+</sup> Ap ln Re<sup>∗</sup> where Re<sup>∗</sup> <sup>=</sup> u∗δ ν and AP can be theoretically expressed as AP = A1[(2p − 1)!!] 1/p, as per Gaussian statistics.

Figure 7 shows the variance in streamwise velocity comparing the inner and outer bends of a sinuous river channel. Both the profiles tend to follow a logarithmic profile after a depth of z<sup>+</sup> > 1000 and tend to approach zero asymptotically after a particular depth (z<sup>+</sup> > 4000). The magnitude of outer bend variance is higher when compared to the inner bend of the channel. This point indicates the erosional behaviour in a sinuous bend. The outer bend of the sinuous channel possesses more magnitude in fluctuating velocity from the mean velocity than the inner bend. However, the mean velocity of the inner bend is more in magnitude than the outer bend, as shown in Figures 3 and 4. This point indicates that outer bend erosion mostly depends on the velocity fluctuation or deviation from mean velocity rather than the mean velocity. The constants for Equation (7) are: (A1, B1) = (3.875, 36.61) for outer bend and (A1, B1) = (1.74, 16.024) for inner bend of the channel. From this, we can understand that as we move from outer to the inner bend, the constants also decrease. This point indicates that the logarithmic law constants are dependent on the flow impact location of the sinuous channel.

**Figure 7.** Streamwise velocity variance in turbulent boundary layers for flow at inner and outer bends of the sinuous river channel.

The relationship between structure function and outer bend erosion can be further established by analyzing the higher-order velocity moments (Equation (8)). Before investigating the higher-order moments of velocity fluctuations, we have to check the convergence of higher-order moments. It was achieved by multiplying the marginal probability density function (PDF) by the velocity fluctuation moments in the near-bed region. Marginal PDF is defined as the probability of values of continuous random variable (say P) without referring to the values of the other variable (say Q). In the present study, marginalized PDF of streamwise velocity (u) is considered. Figure 8a,b show premultiplied PDFs for 2p = 2 and 4, respectively at flow depth z/h = 0.08. Here, we notice that the area covered by the moments of order 2p = 2 and 2p = 4 for both outer and inner bend, respectively, are captured by the available data. In other words, we can say that there is convergence at higher-order moments. Furthermore, the premultiplied velocity fluctuation PDFs is greater in the case of outer bend, indicating that the outer bend suffers from erosion.

Higher-order moments for 2p = 4 and 6 were calculated with respect to flow depth. Figure 8c shows higher-order moments for 2p = 4. The result is similar to that of 2p = 2, where both outer and inner bend profiles follow logarithmic law. The outer bend fluctuation is more as compared to inner bend, which indicates more erosion chances in the outer bend. The fitted constants values for 2p = 4 are (A2, B2) = (8.569, 80.962) for the outer bend and (A2, B2) = (3.986, 35.091) for the inner bend of the sinuous channel. This result was similar to 2p = 2 moments where constants are greater in the outer bend of the channel. The magnitude of higher-order moments was also compared for outer and inner bends of the sinuous channel. Figure 8d,e show higher-order moments for 2p = 2, 4, and 6 at different flow depths for outer and inner bends in a sinuous channel, respectively. We observe that as the order of moments increases, the magnitude also tends to increase. Both the profiles for outer and inner bends at different order moments tend to follow the logarithmic nature and asymptotically tend to zero after some distance. For higher moment, i.e., 2p = 8, 10, etc., it follows the same trend.

**Figure 8.** Premultiplied probability density functions (PDF) of normalized velocity fluctuations u<sup>+</sup>2pP(u+) at z/h = 0.08 with moments (**a**) 2p = 2 and (**b**) 2p = 4. (**c**) Moments of order 2p = 4 for streamwise velocity as a function of wall-normal distance. Moments of different orders of streamwise velocity fluctuation as a function of wall normal distance for flow subjected to (**d**) outer bend and (**e**) inner bend.

Further investigation was carried to see the flatness or kurtosis in the inner and outer bend of the channel (Figure 9). The flatness factor of streamwise velocity was calculated as:

$$\mathbf{F}\_4 = \frac{\left< \mathbf{u}^{+4} \right>}{<\mathbf{u}^{+2}>^2} \tag{9}$$

**Figure 9.** Flatness factor as a function of the wall distance for outer and inner bend of the sinuous channel.

The flatness factor in the outer bend of the channel is greater than that in the inner bend. From Figure 9, we can comment that from <sup>z</sup><sup>+</sup> <sup>≈</sup> 1000, the flatness in the outer bend has increased in comparison to the inner bend. Most points in the outer bend follow F4 > 3 which represents distribution with a peaky signal characteristic. On the other hand, points in the inner bend follow F4 < 3 representing distribution with a flat characteristic. The average kurtosis of the outer bend is 3.06, and that of inner bend is 2.91 throughout the flow depth. This characteristic indicates the erosional behaviour in the outer bend.

## **4. Conclusions**

Higher-order moments of streamwise velocity fluctuations were studied to explain the erosional behaviour of the outer bend of a sinuous river channel. The higher-order moments in association with the morphological changes of a channel contribute a novel approach to understand the turbulent nature in a sinuous channel. The stresses are found higher in the outer bend, which are well reported earlier in the literature. Analyzing the high-order velocity fluctuations has provided clear insights into the scour mechanism near the outer wall of the sinuous bend. The profile of the velocity fluctuations for inner and outer bends at higher moments 2p = 2, 4, and 6 show logarithmic nature from a particular depth (z<sup>+</sup> > 1000). In all these cases, the magnitude of outer bend velocity fluctuation is more when compared to the inner bend. Premultiplied PDFs for 2p = 2 and 4 for outer bend are greater than those of the inner bend. Both the points indicate that though the mean velocity at inner bend is higher than the outer bend, its erosional behaviour mainly depends on the fluctuations from the mean velocity. Further, the constants Ap, Bp of the logarithmic law for velocity fluctuations depend on the location of the flow impact in the sinuous channel. The flatness factor or kurtosis of streamwise velocity was also found to be slightly higher in the outer bend as compared to inner bend. A similar investigation can be carried out with field data. The analysis of structure function can also be applicable to different sediment particles, flow condition, and numerical modeling.

**Author Contributions:** J.T. and J.B. did the experimentation, analysis, and wrote the first draft of the manuscript; B.K. supervised the work; J.T., J.B., B.K., and G.O. completed the final draft. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

**Antonija Harasti \* , Gordon Gilja , Kristina Potoˇcki and Martina Lacko**

Department of Hydroscience and Engineering, Faculty of Civil Engineering, University of Zagreb, 10000 Zagreb, Croatia; gordon.gilja@grad.unizg.hr (G.G.); kristina.potocki@grad.unizg.hr (K.P.); martina.lacko@grad.unizg.hr (M.L.)

**\*** Correspondence: antonija.harasti@grad.unizg.hr

**Abstract:** Bridge piers on large rivers are often protected from scouring using launchable stone, such as a riprap sloping structure. While such scour countermeasures are effective for pier protection, they significantly alter flow conditions in the bridge opening by overtopping flow and flow contraction, deflecting the formation of the scour hole downstream and exposing the downstream riverbed to additional scour. This paper provides a comprehensive and relevant review of bridge scour estimation methods for piers with a riprap sloping structure installed as a scour countermeasure. Research on empirical methods for bridge scour estimation is reviewed and analyzed with formulae used for comparable structures—complex pier formulae and formulae for river training structures. A summary of relevant formulae applicable to piers with installed scour countermeasures is provided, as well as a discussion on the possible future research directions that could contribute to the field.

**Keywords:** bridge scour; empirical formulae; riprap sloping structure; flow contraction; overtopping flow

## **1. Introduction**

The majority of bridges have been built to provide an effective connection between the banks over the waterways, impacting society both economically and politically, and at the same time interacting with the waterway flow regime [1]. During the lifespan of the bridge crossing rivers, changes in the flow and sediment regime are anticipated, as well as consequent change in the variable action-imposed loads on the structure [2], or even complete undermining of the foundation soil [3], which can result in bridge failure or even collapse. Economic losses resulting from traffic disruption following the bridge failure exceed its construction value [4], making bridges critical infrastructure assets [5].

The bridge failures data recorded world-wide indicate that hydraulic causes (scour, floods, stream instability, lateral migration and floating debris) are prevailing among the factors causing bridge collapses: Imhof [6] reports that flooding/scour was the most frequent natural hazard (66%) causing bridge failure in Europe and North America; Muñoz Diaz et al. reported scour as causing 35% of bridge failures in Columbia, with an additional 7% failing from overtopping/floods; Schaap and Caner [7] reported 45% of failures in Turkey were hydraulic-caused, and 22% were directly caused by scour; available data for United States suggests that hydraulic-caused failures account for >50% of the collapses, and scour >20% [8–11]. Common materials used for construction of bridges are concrete/reinforced concrete or steel, with a 75-year design life in service [12]. According to the US bridge failure data, the majority of failed bridges are steel bridges (>60%), although their portion in the National Bridge Inventory is significantly smaller than the concrete ones, at 30% compared with 65%, respectively [13]. The average age of the bridges at the time of the failure was 64 years, and <50 years for steel and concrete material types, respectively.

Among the various types of erosion processes that occur in the riverine environment, the ones influencing bridges are generally divided into long-term general scour, contraction

**Citation:** Harasti, A.; Gilja, G.; Potoˇcki, K.; Lacko, M. Scour at Bridge Piers Protected by the Riprap Sloping Structure: A Review. *Water* **2021**, *13*, 3606. https://doi.org/10.3390/ w13243606

Academic Editor: Mouldi Ben Meftah

Received: 17 November 2021 Accepted: 13 December 2021 Published: 15 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

scour and local scour [14], acting independently or in combination with other hydraulic causes. Bridge piers obstruct the flow, forcing it to accelerate around the pier and generate large-scale turbulence structures, consequently increasing the turbulence in the flow downward towards the bed (horseshoe vortex) and flow downstream behind the pier (wake vortex) [15], causing local scouring. Therefore, local scour is inherently associated with the hydraulic structures interfering with the natural flow field [16] and its occurrence is often simultaneous with on one (or both) of other two types of scour, causing deepening of the entire riverbed in the vicinity of piers and abutments.

Bridge piers can be protected from scouring by either indirect or direct methods [17]. Indirect methods, such as collars [18], vanes [19], sacrificial piers [20], or self-protected piers with openings within the pier [21], change the flow pattern around the bridge pier and thus reduce the strength of downflow and horseshoe vortex [22]. Direct methods (riprap, geobags, cable-tied blocks or similar launchable material) increase the riverbed resistance in order to withstand the turbulent flow around the piers [23]. They are less expensive and easier to construct in comparison to indirect methods [24], and flexible which allows them to be retrofitted if necessary. The most common bed-armoring protection against scour is riprap [25], which consists of a material larger and heavier than riverbed sediment with the main purpose of preventing removal of the riverbed sediment downstream [26]. There are several types of riprap countermeasures: horizontal riprap layer [17] placed just above the riverbed level, riprap sloping structure [27] formed by stones conically mounded around the pier, gabion mattresses [28], sacks filled with stones [29], etc. Although riprap is considered stable, there is evidence of riverbed particle erosion during floods and consequent riprap slumping and sliding. Breusers et al. [30] pointed out that riprap protection can also induce scour, thus not fully fulfilling its purpose. Nielsen et al. [31] investigated the velocity distribution and critical bed shear stress under the riprap protection. Fredsoe et al. [25] described the horseshoe vortex as the primary mechanism that causes undermining the stones at the junction between the riprap and the riverbed. This way, riprap protection causes formation of the deflected scour hole next to the riprap scour protection, introducing hazard to the pier or adjacent river structures [32].

Scouring locally deepens the riverbed next to bridge piers, exposing a larger area of the pier to direct oncoming flow, thus subjecting them to greater hydrodynamic forces in comparison to nonscoured conditions [33]. Changes in hydrodynamic loading coinciding with other natural hazards, e.g., earthquakes, can have devastating consequences for both infrastructure and society [34]. Morphodynamic analysis of the river reach influencing the bridge is integrated into the bridge management system [35] in order to log changes developed since the historical reports and assess the current condition of the riverbed. Riverbed condition assessment is not unified—each county has adapted inspection procedures to the diversity of their infrastructure, respecting viable methods and instruments for gathering relevant and reliable data [36].

Scour holes reach their maximum depth during floods [37], after which they are backfilled by the bedload of the receding flow, reducing the apparent scour depth [33,38–40]. Most research in recent decades has focused on describing scour depth based on maximum equilibrium scour depth and for steady flow conditions [41,42]. Morphodynamic changes in the riverbed and consequently the scour process around bridge structures, when considered on a larger spatial scale, are not only caused by anthropogenic changes in the watershed, but also by changes in the natural flood regime and extreme climatic events such as floods, that affect the overall stability of bridge structures. Therefore, bathymetric data acquired post-flood are not reliable enough to determine the full extent of scouring, and analyses based on such data can lead to erroneous estimates of bridge safety. Realtime monitoring of scour development during floods is recognized as essential [43] to enable safety monitoring and systematic inspection of the riverbed condition and avoid unnecessary maintenance. Some of the monitoring equipment used for scour detection are magnetic sliding collars [44], sonar devices [45], float-out devices [46] or tilt sensors [47]. The current scour-monitoring systems are being constantly enhanced through development of advanced monitoring techniques, mainly to achieve more robust devices, reduce maintenance cost and enable integration with other sensors. Recent advances in scour monitoring include assessment based on vehicle-induced vibrations measured by acceleration sensors [48], influence of soil characteristics on the change in the predominant natural frequency of a bridge pier [49], scour depth measurement using vibration-based microelectro-mechanical systems [50], instrumented particles [51], interferometric synthetic aperture radar stacking techniques [52], amongst others.

Regardless of the data acquirement method, in order to complete the risk assessment for preliminary hazard analysis in the bridge management system, hydraulic analysis must be performed. Hydraulic analysis will combine the scour data with other relevant variables: flow environment, pier geometry and riverbed composition, to quantify the scour hazard in comparison to foundation depth. Flow velocity, flow constriction and scour protection measures are the most important factors that should be considered when assessing scour risk to bridges [53]. For this purpose, empirical formulae are mostly used that calculate the maximum scour depth, taking into account all relevant variables as input. These empirical formulae are reliable if their usage adheres to conditions under which they were developed and validated on field-measured data, if possible.

Previous studies have focused on the scour development for a wide range of pier geometries, arrangements and spacing, riverbed compositions, sediment transport conditions, etc. However, there exists an evident gap in knowledge regarding the scour occurring at piers protected with riprap sloping structure and associated scour at the periphery of the riprap mound. The aim of this paper is to provide a comprehensive and relevant review of bridge scour estimation methods and provide discussion of their applicability for piers with riprap sloping structure installed as a scour countermeasure. Research on empirical methods for bridge scour estimation is reviewed and compared to complex pier formulae and formulae for river training structures similar to riprap sloping structure scour protection. Scour formulae in this paper are categorized considering the following: the data used for their development (field vs. experimental); pier geometry (single vs. complex piers); scour type (local vs. contraction), and riprap type (riprap sloping structure vs. layer riprap). A summary of relevant formulae applicable for piers with installed scour countermeasures is provided, as well as a discussion on the possible future research directions.

## **2. Bridge Scour Assessment**

Most often, the approach for developing a scour depth formula is the application of dimensional analysis technique to an experimental dataset. The overall formula is assessed by validation with field-measured data and incorporation of one of the mathematical fitting models. Dataset for establishing formula can be obtained in the field and experimental environment.

## *2.1. Field and Experimental Data*

Compared to the field data, the experimental data obtained in the laboratory provide a controlled flow environment, which is why the maturity of the scour hole is reliable only from the laboratory experiments [54]. However, the experimentally developed formulae have numerous disadvantages: a limited applicability in a range of input data due to scaling issues between geometric and sediment scales; oversimplified conditions that do not correspond to the field conditions (rectangular flume, uniform sediment gradation, clear-water conditions, and steady uniform flow), and time-consuming simulations to acquire the desired dataset. Time to reach equilibrium scour depth condition is reached asymptotically, so duration of the experimental test needs to be much more than at least 48 h to reach at least 90% of the equilibrium [55,56]. Melville and Chiew [43] define a formula for time needed for development of the equilibrium scour hole as s rate of 5% of the pier diameter in a 24 h period, which means that based on first 24 h-long experiment duration, equilibrium scour depth can be calculated. Their results showed that after 10% of the total time, scour depth is in the range of 50% to 80% of the equilibrium scour depth. To

shorten time-consuming simulations, some authors accept semiequilibrium conditions, so each experiment lasted 4 to 6 h, and subsequently extend the scour depth values to 100% of the equilibrium scour depths.

Sheppard and Melville [55] noticed why scour depths are overestimated in small-scale experiments—because the forces of the pressure gradients on sediment particles near the smaller bridge piers are much greater than for larger piers. A similar observation was made by Huang et al. [57] in a scale effects study, where a comparison of results from a large-scale numerical model and a small-scale physical model was conducted. They found better prediction of turbulence flow pattern and sediment scour by the numerical model in comparison to the physical model, because the physical model introduces Froude similarity and ignores the effect of turbulent Reynolds number.

The numerical model or computational fluid dynamic (CFD) model is powerful tool to simulate interaction between the sediment transport process and the vortex system in the complex geometry environment. In contrast to time-consuming laboratory experiments and field investigations, numerical simulations are inexpensive and fast way to collect a large and comprehensive amount of the essential scour data. Numerous turbulent models have been already developed to calculate turbulent velocity and sediment scour (*k* − *ε*, *k* − *ω*, LES, RSM, RNG and others). For example, Aly and Dougherty [58] examined different turbulence models in predicting the bed shear stress around different pier scour countermeasures. They concluded that *k* − *ε* performed better than the *k* − *ω* model in approximating the experimental data, since *k* − *ω* simulates adverse pressure gradients, while the *k* − *ε* realizable model successfully simulates velocity components. Alemi et al. [59] used the LES turbulence model to investigate the flow around a complex bridge pier on the scoured fixed bed. They validated the numerical model based on laboratory experiments by Graf and Istiarto [60] and achieved high correspondence with their results. By combining LES and wall function, a fast simulation process is achieved since the wall function simplifies fully developed flow with zero velocity near the wall region. Developing a scour numerical model can be challenging due to establishing mesh on complex riverbed geometry, which changes gradually as scour occurs. Once the computational grid is established, the flow field will scour the bed form from the first timestep, and the computational domain needs to be re-meshed. Zhu and Liu [61] presented three re-meshing strategies: *z*-level grid, shaved grid, and *σ* grid. Jalal and Hassan [62] examined various types and sizes of cells to identify optimum cell (5–10 mm) that balances the accuracy of results and reduction in computation time. Song et al. [63] developed a 3D scour model, ibScourFoam, based on an immersed boundary method, validated against available data from the literature and their own flume experiment. Their model showed advantage over previous scour models in accurate simulation of the scour process around complex structures.

Experimental and numerical data must be validated by field data to ensure the reliability of experimental results [64]. Field measurements of hydraulic parameters during a flood are essential to understand the morphodynamic evolution of the local scour phenomenon. Despite today's sophisticated survey techniques, it is not easy to accurately measure velocity and depth during floods in real-time due to instrument functionality under oscillating flow depths and instrument durability, as it is exposed to impact damage by floating debris. Scour monitoring is challenging from the perspective of instrument placement, where approximate location of the scour hole needs to be known. Surveying implies collecting data at a certain point in time, which means that the maturity of the local scour hole is unknown at the time of measurement. In addition to the difficulty of obtaining accurate scour depth data, the flow field upstream and downstream of the bridge pier must also be collected. Maneuvering a boat during a flood is difficult and hazardous in vicinity of the structures. Moreover, there are some difficulties in estimation of bed material size which used to be highly variable on site, and inaccuracy due to neglecting scour-hole backfilling effect [65].

For these reasons, in the present state-of-the-art data from field surveys are limited and do not encompass a sufficient range of hydraulic properties and characteristics of the river

channel to represent a solid basis for high-quality research. Therefore, most of the local scour formulae have been developed based on laboratory experiments [14,22,30,66,67]. Yet, there are several empirical formulae developed from field data only. Governing parameters for scour depth (*ds*) formula developed from field data are width of bridge pier (*b*), water depth (*y*0) and Froude number (*Fr*) at the approach section and correction factor for pier-nose shape (*K*1) as presented in Table 1.


**Table 1.** Local scour data obtained by field survey and corresponding formulae.

Many authors have expressed concern about the lack of information available in the published studies regarding the source of dataset used to develop the empirical formulae [55,59,71] and insufficient information about the conditions that were present during the measurements. For example, there is no information on the methods that Zhuravlyov [72] used to obtain the data, neither the dimensions of the structures nor their shapes, and a field survey by Froehlich [68] does not contain a sediment gradation parameter.

## *2.2. Variables Governing the Scour Process*

Numerous studies have aimed to evaluate the accuracy of the relevant formulae for scour depth estimation [73–77]. These comparative studies generally summarize that most existing formulae tend to overestimate the scour depth equilibrium. However, evaluating reliability of scour formulae must be done with caution because single variables within each scour depth formula may be defined with several curves, which implies dependency of a third variable [73].

Scour can be greatly affected by sediment size distribution which can be uniform or nonuniform. Initial movement of smaller particle sizes could be achieved at lower velocities, so scour may occur more readily in sandy riverbeds than in coarser sediment material [58]. For nonuniform sediment environments bed armoring occurs where the coarser gravel size protects finer particles from being transported by the flow. Pandey et al. [78] emphasize that the influence of the armored layer regarding local scour is poorly investigated because sediment transport in the armored riverbed is generally challenging to explicitly define. However, the authors analyzed how the formation of the armored layer in the scour hole stops further scouring processes during clear-water conditions. They proposed a new graphical approach for estimating maximum scour depth in an armor riverbed, showing that changes in dimensionless scour depth (*ds*/ *<sup>b</sup>*2·*y*<sup>0</sup> 1/3) with armor ratio (*d*50*a*/*d*50) are minor when it exceeds 0.5, where *d*50*<sup>a</sup>* is the median diameter of armor layer *d*<sup>50</sup> is median sediment size.

Most of scour formulae differentiate clear-water or live-bed conditions respective to bedload sediment transport: under clear-water conditions no bedload transport is present, while live-bed condition considers initiated sediment bedload transport. Landers and Mueller [73] evaluated the scour formulae by Gao et al. [69], comparing them to field measurements, and concluded that there was underprediction in deeper live-bed scour conditions (for *d*<sup>50</sup> > 3 mm). On the contrary, the HEC-18 [41] formula overestimates scour depth under live-bed conditions, especially for greater flow intensities, while it shows good agreement with clear-water threshold scour [79]. Scour depth at the clear-water scour is considered to be 10% greater than the scour depth in live-bed conditions [58], but

that finding is based on small-scale models. However, recent studies have shown that the reverse phenomena may occur in the field or in large-scale numerical models [79]. Sheppard and Melville [55] assembled field and experimental datasets of scour data from different authors in the literature and eliminated unreasonable ones by quality control screening, where 441 are experimental and 791 are field data. These data were used to test the accuracy of different scour formulae. The best performance showed the new Sheppard– Melville formula which was originated from modification of Sheppard and Miller's [67] and Melville's formulae [80]. Hereafter, Yang et al. [81] proposed a modified Sheppard– Melville formula regarding the influence of flow shallowness ratio (*y*0/*D*), where *D* is pier diameter. They prove the considerable influence of flow shallowness ratio scour on scour depth, which increases with an increasing (*y*0/*D*). Although modified the Sheppard– Melville formula still retains a small level of overestimation, Yang et al. demonstrate better accuracy of the formula than the HEC-18 formula. Overestimation is greater with smaller shallowness (*y*0/*D*) ratio. Yang et al. collected live-bed field data during short-peak flood waves and concluded that live-bed scour peak may exceed the threshold peak.

Changes in the natural flow regime of a river can be defined in terms of ecological [82] and morphological processes [83], by magnitude, frequency, timing, duration, and unsteadiness. In contrast to engineering design practice, where typical flood events with a return period of 100 or 200 years are chosen to account for the scour susceptibility of bridges [84], Tubaldi et al. [85] propose a more rigorous assessment of scour susceptibility, that examines not only a single flood event but also subsequent events, as floods of short duration or low intensity can lead to partial erosion of the riverbed, making it easier for a subsequent flood to reach a maximum scour depth. Therefore, although the flow magnitude is the main factor influencing scour depth development, hydrograph duration and shape also have an important influence on the scour process and should be taken into account. The maximum scour depth is significantly lower when actual hydrograph shape is considered than the equilibrium scour depth for the constant flow rate. According to Plumb et al. [83], the effects of hydrograph shape, addressing the number of cycled hydrographs and duration of each hydrograph, as well as changes in the flow regime, should be examined in the scour development analysis. In addition, the inclusion of more complex sediment grain size distributions in the estimation of scour next to the bridges, along with the hydrograph shape parameters, results in a potentially higher threshold for bed motion state, complicating the prediction of instantaneous flux using formulae for steady-flow conditions [86,87]. Some empirical formulae have been developed for evaluation of time-dependent scour depth, i.e., considering hydrograph shape characteristics [43,88]. However, these formulae were mostly developed for steady-state flow conditions or adapted to unsteady conditions using various approaches, such as superposition of the hydrograph as a sequence of steady-state discharge steps [89–91], introducing a mathematical function for form of the hydrograph [92], or defining a dimensionless effective flow [93,94]. Apart from the fact that real-time field measurements are rarely available at high turbidity during a flood, the main challenge in studying the effects of flow characteristics, i.e., unsteady flow characteristics and the influence of multiple flood waves [95] on bedload transport and the scour process, is to separate the effects of flow conditions from the effects of bed material properties and sediment supply conditions [96]. It can be concluded that in future research it is necessary to investigate the effects on the development of maximum scour depth not only during a single flood event and the associated flood wave characteristics, but also previous flood events should be considered. Particular attention should be paid to the development of formulae and field measurement techniques that would allow the determination and measurement of the temporal evolution of scour depth.

## *2.3. Local Scour Formulae for Complex Piers*

Many studies of local scour formulae refer to the single pier formula [30,70,97] which is a function of pier width (*b*), flow depth (*y*0), Froude number (*Fr*), critical flow velocity (*vc*), and sediment bed size (*d*50). Sometimes the formulae are expanded to include the

shape pier factor (*K*1), the angle of flow attack (*α*) (if the pier is skewed), or the sediment characteristic (*ρs*) [98,99] and generally the formula can be written as follows:

$$d\_s = f\left[\text{flow}(y, Fr, v\_c); \text{ pier characteristic}(b, l, a); \text{ bed material}(d\_{50}, \rho\_s)\right]. \tag{5}$$

Physical, economic, and geotechnical considerations usually indicate the need for bridge piers to be constructed with a column founded on a pile group with a pile cap or on a caisson, so-called complex piers [100]. The process of scour development at complex piers differs from that on uniform piers because lee-wake vortices dominate sediment transport, whereas on uniform piers horseshoe vortex has governing influence on scouring [101]. In recent decades, a number of studies on scouring around complex bridge piers have been conducted. Melville and Raudkivi [102] proposed an approach to calculate the scour around complex piers composed of a cylindrical pier (with diameter *D*) and of a foundation (with diameter *D*∗ and with the top elevation *H*). In their study, an alternative approach using the effective diameter (*De*) is employed. Effective diameter represents a diameter of a circular pile that induces the same scour depth as the scour depth at the complex pier. Melville and Raudkivi [102] examined the influence of effective diameter concept and proposed following formula:

$$D\_{\varepsilon} = \frac{D \cdot (y\_0 + H) + D \ast \cdot (d\_s - H)}{d\_s + y\_0} \tag{6}$$

Finally, Melville and Raudkivi's scour depth formula for complex piers introduces three scour zones dependent of the top of the foundation elevation, as presented in Table 2.

**Table 2.** Scour depth formula around complex piers based on influence of the top of the foundation.


<sup>1</sup> Except *<sup>H</sup>* < 0.7·*<sup>D</sup>* and *<sup>D</sup>*/*D*<sup>∗</sup> < 0.6

Even though using the equivalent single pier simplification shortens the computational time considering the complex flow field and turbulence around the pier, Melville and Raudkivi [102] showed that it could lead to conservative estimates of scour depth.

One of the most commonly used [65,75,103,104] scour formulae is the HEC-18 formula, developed by the Federal Highway Administration (FHWA) [41]. To account for the effects of foundation geometry on scour depth, the superposition of scour components is proposed. In this method, each pier component is calculated separately and superimposed on the total scour depth, as follows:

$$d\_s = d\_{s,pier} + d\_{s,\text{ pile }cap} + d\_{s,\text{pile }ground} \tag{10}$$

Afterwards, Coleman [105] expanded scour estimation of complex piers (comprising column, pile cap and pile group) into the five following cases dependent of the top of the foundation elevation: (I) buried pile cap; (II) pile cap at the bed level of the scour hole with no exposure of the pile group; (III) pile cap and pile group extending over the base of the scour hole while water level is above the pile cap; (IV) pile cap and pile group extending over the base of the scour hole while water level is at the top of the pile cap; and (V) pile cap and pile group extending over the base of the scour hole while water level is below the pile cap. Coleman also provides formulae for equivalent pier diameter (*be*) based on variation of the pile cap elevation (*H*) relative to the bed level (Table 3.):


**Table 3.** Equivalent pier diameter formulae based on influence of the pile cap elevation.

A similar study by Sheppard and Glasser [66] resulted in a detailed iterative procedure for determining effective diameter for the geometry case where the top of the pile cap is above the bed level with no exposure of the pile group—equivalent to the piers founded on the exposed caisson. The outlined procedure was developed in clear-water and livebed conditions, whereas effective diameter is a function of the shapes, locations, and orientations of each pier component, and as a function of the sediment and flow properties (Table 4).

**Table 4.** Procedure for calculating effective width of complex piers for estimating scour depth.


<sup>1</sup> Caisson is interpreted as buried pile cap.

Furthermore, Jannaty et al. [77] conducted an investigation to determine the cause of the Adinan complex bridge failure. They examined the performance of several formulae and demonstrated substantial distinction between calculated and measured scour depth values. Several stages of scouring can be noticed in their study: (I) if the bottom of scour hole is above the up edge of the pile cap, the pile cap has a protective effect against scour; (II) if the bottom of scour hole is between the up and the down edge of the pile cap, the column is the prominent component to produce scour; (III) if the scour hole induced by the column reaches the down edge of the pile cap, the pile cap will play the most important role in scouring. Analysis of results show that the cause of Adinan Bridge failure is the large width of the foundation in comparison to the column width. Additionally, scouring was intensified due to lot of accumulated debris. Based on study by Jannaty et al., it can be concluded that scour depth increases if obstruction area (either foundation or debris) increases. This leads us to the conclusion that in estimating scour depth, the component of flow area reduction should be considered in addition to local scour. Similarly, Yang

et al. [106] found that for complex bridge piers in close proximity, the upstream pier significantly affects scouring at the downstream pier in comparison to individual piers, reducing scour rate for clear-water conditions and increasing the bedform celerity for live-bed conditions.

## *2.4. Contraction Scour Formulae*

One of the first contraction scour formulae was presented by Laursen and Toch [107]. They evolve Straub's contraction scour for the long contractions to the local scour at the pier. The formula is valid for bridges where the spacing between the piers is much larger than the piers' width (for contraction of cross-section approximately 10%). The effect of contraction scour will become present when scour depth achieves a value that corresponds to scour depth which would occur in a long contraction, or when scour holes from adjacent piers overlap. The formula is as follows:

$$\frac{d\_{s,\mathcal{L}}}{y\_0} = \frac{1}{\left(1-\beta\right)^{9/14}} - 1,\tag{18}$$

where *β* is width ratio of the contracted to the uncontracted profile. Briaud et al. [108] verified contraction scour data using previous experimental results, employed a new parameter *v*1, that presents averaged velocity at the contracted section, and provided the following formula:

$$\frac{d\_{s,c}}{y\_0} = 2.21 \cdot \left( 1.31 \frac{v\_1}{\sqrt{\mathcal{G}\mathcal{Y}\_0}} - \frac{v\_c}{\sqrt{\mathcal{G}\mathcal{Y}\_0}} \right). \tag{19}$$

Furthermore, one of the most commonly used contraction scour formulae is HEC-18 [109], that assumes contraction scour as a long contraction, where the length of the contracted section is longer than the section width. The HEC-18 contraction scour formula is developed as a part of total scour formula, which presents superposition of both separately calculated processes. The formula contains flow, bottom channel width, and Manning's roughness coefficient ratios as provided below:

$$\frac{d\_{s,c}}{y\_0} = \left(\frac{Q\_1}{Q\_0}\right)^{\frac{6}{7}} \left(\frac{L\_0}{L\_1}\right)^{k\_1} \left(\frac{n\_1}{n\_0}\right)^{k\_2} - 1. \tag{20}$$

Like many other contraction scour formulae in the literature, HEC-18 also overestimates the scour depth [75,110]. One of the possible reasons for that inaccuracy could be neglecting the possible dependence of local and contraction scour [65].

Recently, several researchers have observed uncertainty in the scour formulae due to the interaction of contraction and local pier scour estimation. Hong [110] conducted laboratory experiments comparing the bridge cross-section with piers and abutments with the cross-section without piers. The results demonstrate that the model without piers has 25% greater scour depth. A decrease in contraction scour depth implies possible interactions between local scour and contraction scour processes. Contraction scour development requires a longer time to reach equilibrium than local scour, as well as a possible consequence of local and contraction scour dependance. Additionally, experimental results of contraction scour depth were compared to the values obtained by Laursen's formula and results shows overestimation of Laursen's formula results of about 30–60%. Recently, Saha et al. [104] carried out laboratory experiments in the presence of flow contraction and compared results with reference local scour depth data obtained by new theoretical method developed from a combination of HEC-18 and M-S formulae. Additional contraction scour estimation showed that the discharge contraction ratio *q*2/*q*<sup>1</sup> has a great influence on the contraction effect on the total scour depth, where the cross-section at the upstream face of bridge is labeled with number 2 and approach section with number 1. However, their experiments were only conducted under clear-water conditions, and the nonuniform sediment size is neglected. Mueller and Wagner [64] emphasize that long-term in situ measurements are

essential to separate local scour from contraction to identify reference bed elevation. They suggested that long-term scour in the uncontracted section can be determined by adding the general slope of the line to eliminate general aggradation or degradation. Afterwards, short-term scour determination should be performed by comparing the uncontracted and contracted cross-sections in the preflood and flood conditions.

## **3. Riprap Sloping Structure**

When riprap is selected as a scour countermeasure, it is usually placed flushed with the riverbed level if flow conditions are favorable for detailed installation, or as a riprap sloping structure when flow conditions do not allow fine maneuvering of the machinery.

Some authors consider layer riprap as a more suitable protection than a riprap sloping structure [26,30,107] because of excessive exposure of the sloping structure that induces unwanted contraction. Recommended placement of the layer riprap depends on the hydraulic environment and sediment transport, i.e., the occurrence of bedforms. Riprap should be placed deep enough so it does not protrude above the bed level and disturb the flow. Installation in rivers with significant amounts of bedload should be buried below the estimated through of the dunes occurring during floods [10]. The number of studies in the literature describing layer riprap and its failure mechanisms is higher than that of riprap sloping structures, which have been rarely investigated. Generally, hydraulic processes that occur downstream of the riprap sloping structure and its consequences are ambiguous—flow obstruction, debris accumulation and scouring at the periphery of the riprap mound. Since there is no such research about estimating scouring around the riprap sloping structure, an analogy with complex piers can be drawn. Contours of the piers protected with riprap sloping structure (pier + riprap) are similar to the contours of complex piers (pier + caisson), as can be seen in Figure 1.

**Figure 1.** Similarity of outer contours between the riprap sloping structure (**left**) and complex piers (**right**), outlined in red.

If analogy between the riprap sloping structure shape and complex piers geometry is considered valid, then the scour formulae applicable for the riprap sloping structure also have to be case-specific depending on the top elevation of the riprap, similar to the scour formulae for complex piers that are distinguished depending on the caisson submergence [66,77,102,105]. According to Melville [111], there are three cases for scour around complex piers: case I—where the top of the caisson foundation is below general bed level; case II—where the top of the caisson foundation is above general bed level; case III—where the top of the caisson foundation is above the water surface, as depicted by Figure 2.

**Figure 2.** Pier scour countermeasure retrofitting cases distinctive by the position of the top of the riprap sloping structure relative to the riverbed level: top of the riprap below the general bed level (**left**); top of the riprap above the general bed level (**middle**); and top of the riprap above the mean water level (**right**).

Melville [111] and Ghorbani [112] draw further conclusions: if the top of the caisson (or riprap sloping structure) is exposed above general bed level, the scour depth is increased; if the overflow above the top of the pile cap (crest of the riprap) is less shallow, the scour depth is increased. The riprap sloping structure is a wide structure, in terms of reducing the flow area, so conclusions about wide piers might be comparable. Yang et al. [81] stated that the downflow in front of the wider piers (*y*0/*D* ≤ 1.4) is weakened, as well as there being horseshoe vortex. Thus, wake vortices consequently remain the main turbulence structures that produce scour hole, the maximum of which will occur at side of the pier.

The design of riprap considers calculating the size of riprap stone that can withstand flow attack. One of the first formulae for riprap sizing was developed by Izbash [113]:

$$d\_{I} = y\_{0} \cdot 0.347 \cdot \frac{Fr^{3}}{(S\_{I} - 1)}.\tag{21}$$

As can be seen from previous formula, water depth (*y*0), Froude number (*Fr*), and specific gravity of riprap (*Sr*) were governing parameters in riprap sizing formulae [30,114,115], until Parola [116] introduced new characteristic factors. Factors are proportional to the size of the riprap stones (*Kr*) and related to pier geometry—pier shape (*Ks*) and pier width (*Kb*). Subsequently, several formulae with pier characteristic factors were developed [46,117]. Furthermore, Froehlich [118] compared the range of riprap sizes obtained by formulae from different authors, and came out with conclusion of quite a wide range of results, and therefore proposed a new formula for calculating minimum diameter of loose rock riprap, introducing new parameters that introduce crossflow shear (*Kw*), transverse pier spacing factor (*Kp*) and approach flow influence (*Kα*) as follows:

$$d\_r = \mathcal{y} \cdot \mathcal{K}\_r \cdot \mathcal{K}\_b \cdot \mathcal{K}\_{\omega} \cdot \mathcal{K}\_p \cdot \mathcal{K}\_s \cdot \mathcal{K}\_a \cdot \mathcal{F} r^3. \tag{22}$$

Besides riprap sizing, the design of the riprap sloping structure stability includes calculating an appropriate slope to prevent the riprap stones from sliding into the scour hole, considering protection at the toe due to undermining, and determining the riprap stone diameter to resist the hydrodynamic force of the flow. For riprap side slope, Park et al. [27] used the U.S. Army Corps of Engineering (USACE) method for calculating riprap size as follows:

$$d\_{r, \text{j0\%}} = y\_0 \cdot \frac{1.1 \cdot \mathbb{C}\_{\text{s}} \cdot \mathbb{C}\_{\text{v}} \cdot \mathbb{C}\_{\text{t}}}{K\_{\text{sl}} \cdot 1.25 \cdot (S\_{\text{r}} - 1)^{1.25}} \cdot Fr^{2.5} \,, \tag{23}$$

which was developed primarily for bank and channel protection. Parameter *dr*,30% represents particle diameter, corresponding to a 30% finer grain-size of the granulometric curve, *Cs* is stability coefficient for incipient failure, *Cv* is vertical velocity distribution coefficient, *Ct* is blanket thickness coefficient and *Ksl* is side-slope correction factor. According to Breusers et al. [30], the horizontal dimension of riprap protection depends on the pier diameter, and it needs to be at least two times larger.

## *3.1. Riprap Failure Mechanisms*

Riprap failure mechanisms in clear-water conditions can be divided into: (I) shear failure—dislodging of individual riprap stones due to hydrodynamic forces; (II) edge failure—undermining of the riverbed at the toe of the riprap; and (III) winnowing failure movement of finer material through voids between riprap stones principally initiated by turbulence [117]. In addition, two more failure mechanisms have been identified under live-bed conditions: (IV) bedform-induced failure—fluctuations in flow of the bed level prompted by bed features, such as ripples and dunes, causes lack of stone stability; (V) bed-degradation-induced failure—due to general bed degradation, riprap stones protrude above the general bed level causing a reduction in bed shear stress and consequently disintegration of the riprap structure [119]. Fredsoe et al. [25] described horseshoe vortices in front of the individual stones as the primary reason for undermining at the junction between the riprap countermeasure and the bed in a steady parallel current. The following conclusion is determined—scour depth increases if the Shields parameter increases, a slope of the revetment increases, and if the stone shape is angular. Vasquez et al. [24] conducted physical and numerical investigation of scour around the riprap layer protection in Golden Ears Bridge. In their study, the riprap layer is called riprap apron. The bridge is located in a wide section of a sand-bed river just upstream from a bifurcation, subjected to migrating dunes under live-bed conditions, and with piers oriented at an angle to the flow attack, except the pile group that is exposed above the natural riverbed level. The study proves that a 16 m-wide and 2 m-thick riprap layer successfully eliminates local scour around the pier, even in the presence of an unfavorable combination of passing dunes and general scour. The riprap, however, will partially destabilize and reallocate individual riprap stones, especially at the edges of the layer to the form of a semiconical mound around the pier, but it will still prevent scour hole from forming near piers.

However, at the contacts of the riprap and sediment bed in terms of live-bed conditions, the edges of the riprap will inevitably start to destabilize. The upper stones at the edges roll down the riprap slope or towards the dunes through into the scour hole. Vasquez et al. [24] noticed that after the dune passes by, the riprap side slope formed a conical mound with an angle of repose of about 30◦. This process continues downstream until the riprap becomes scattered and forms a riprap mound. Chiew [120] states that such a riprap mound still serves its function of protecting the pier from erosion, but nevertheless the riprap mound shifts the scour area away from the riprap toe and scour continues downstream. This phenomenon has been noticed in Croatian rivers [121]. Gilja et al. establish a sediment transport model to investigate morphological characteristics of a section of the Sava River where a bridge with three piers remediated with riprap sloping structure are placed. It has been noticed that erosional processes affect downstream sections of the river due to increased flow velocities at the constricted bridge opening. The location of the final scour hole is unknown, because it is formed based on the interaction between the flow and the structure under site-specific conditions [122]. Petersen et al. [123] noticed scour of the riprap in the farther downstream area. They conducted a detailed experimental and field investigation on the three-dimensional flow field via Particle Image Velocimetry (PIV) measurements to improve the understanding of scour in the interaction of the riprap protection and the sand. The study was based on monopiles in the marine environment in terms of steady current, waves and the combination of current and waves, where analogy with sea current and river flow can be correlated. They concluded that a pair of symmetrical, counter-rotating vortices induced by a steady current cause significant scour hole downstream of the scour protection. Whitehouse et al. [71] support the statement that the riprap placed around the pier and above the level of the natural bed initiates the development of the secondary scour response. As the scour continues downstream at the edge of the riprap, a deeper scour hole is created than around the unprotected bridge pier. In addition, it has been shown that with this type of riprap protection, scour wakes extend to distance of 100 times the pier diameter. The riprap scour countermeasure is considered

to be inadequate, because after secondary scour, more vortices appear near the piers and accelerate the scour process [124].

#### *3.2. Similar Structures Analogy*

Groynes, spur dikes, and sloping abutments are geometrically similar to the riprap sloping structures, and thus comparable complex interactions between the flow and riverbed material occur—shallow and on occasion supercritical flow overtopping the structure and hydraulic jump forming on the downstream slope, undermining the toe.

Affecting parameters in local scour formulae at the groyne are same as those for local scour at bridge piers (hydraulic and sediment parameters), except for geometrical parameters that encompass scour depth's dependency on the distance from the first groyne or from the entrance of the bend, as can be seen in Przedwojski [125] formula:

$$\frac{d\_s}{y\_0} = \left[\frac{y\_{\text{un}}}{y\_0} + \frac{Q\_s}{Q\_0}\cos\left(\frac{2\pi}{L\_{\text{m}}}\mathbf{s}\right) + \sin(\mathfrak{a} - \mathfrak{H}0)\right]^{\mu - 2}.\tag{24}$$

Another formula developed by Rashak and Khassaf [126] in clear-water conditions for T-shape submerged groynes is as follows:

$$\frac{d\_s}{y\_0} = 1.489 \cdot e^{(4.117 \cdot \text{Fr})} + 4.117 \cdot e^{(-3.292 \cdot \frac{\text{r}\_0}{\text{yr}})} - 3.292 \cdot e^{(0.002 \cdot \frac{\text{l} \cdot \text{yr}}{\text{yr}})} + 0.002 \cdot e^{0.747 \cdot \frac{\text{l} \cdot \text{r}}{\text{yr}}} + 0.747 \cdot e^{(-1.092 \cdot \text{n}\_3)}.\tag{25}$$

As can be seen in abovementioned formulae, local scour around groynes depend on parameter of spacing between groynes (*sg*). Since one groyne case is comparable to the riprap sloping structure, Rashak and Khassaf's formula can be considered for the riprap sloping structure in the case where spacing between groynes is zero (*sg* = 0). In a sequence of groynes, maximum scour depth occurs at the nose of the upstream groyne.

According to Han and Lin [127], groynes are most vulnerable when the upcoming flood reaches an elevation just above the groyne crest and a shallow submerged condition is reached. Such overtopping flow has a vertical stream direction down the lee side of the groyne, and thus causes small sediment particles to be carried away at relatively low threshold velocities. Rashak and Khassaf [128] noticed that the magnitude of the groyne submergence affects scour depth inversely, which means that an increase in submerged ratio (*y*1/*y*0) will result in smaller scour depth ratios (*ds*/*y*0). During relatively deep submergence, the overtopping flow near the free surface is parallel to the flow direction, while recirculatory motions are present at the bottom around the groyne. Meaning that in the low submergence level case, the recirculation region affects stones at the groyne surface.

McCoy et al. [129] investigated the horseshoe vortex system and shear distribution around the groyne to clarify evolution of the scour hole around the groyne region. They showed that deflection of the horizontal, detached shear layer will be magnified if the groyne crest elevation is higher. Furthermore, they support statement that mass exchange is larger in the submerged conditions than for emerged groynes due to a large increase in the inflow velocities downstream of the groyne. Melville and Coleman [14] mentioned that the vortex system will not change in higher water levels, which means that flow depth has a limited influence on scour depth.

Pandey et al. [130] have tested the accuracy of formulae for scour depth around groynes developed by different authors, concluding that although the position of scour depth varies remarkably, a maximum value typically occurs near the nose of the groyne on the upstream side and spreads up to a width of three times the spur dike length. The extent of the scour region is greater on the downstream side compared to the upstream, but the volume of the upstream scour hole is about 65% of the total volume. They indicate that Froude number is a significant parameter in scour depth around the groyne, and other parameters play secondary role, which was later also confirmed by Rashak and Khassaf [126]. Giglou et al. [131] conducted an analysis of flow pattern around spur dikes by a 3D numerical model, and showed that vortex length behind the spur dike is four

times longer than the length of spur dike, and with approximately 1.2 times the spur dike length. Pandey et al. [132] developed two novel methods to estimate maximum scour depth around the spur dike in a uniform sediment condition that consists of three standalone machine-learning approaches. Experimentally collected data were obtained in a clear-water condition with uniform sediment material. The result of their model is presented as the ratio of scour depth and spur dike (*ds*/*Ls*), while input parameters were ratio of average to critical velocity (*v*0/*vc*), ratio of water depth to spur dike length (*y*0/*Ls*), ratio of spur dike length to mean sediment diameter (*Ls*/*d*50) and densimetric Froude number (*Fd*50). Based on sensitivity analysis of the input parameters, *Ls*/*d*<sup>50</sup> showed the most significant influence on model performance. Eventually, statistical metrics showed the good performance of the developed model in assessing maximum scour depth near the spur dike.

Scour hole geometry at the abutment is similar to those at groynes by means of strong scour activity as a consequence of contraction scouring, which causes a scour hole to be located slightly downstream along the centerline of the watercourse, complementary with turbulent intensity pattern and the time-averaged velocity pattern [109]. However, at long abutments, the scour hole is slightly more elongated to the downstream in comparison to the short abutments [14]. It can be concluded that increasing the obstruction of the flow area will cause larger principal vortex that deflects the scour hole away until it reaches a more elongated shape.

## **4. Discussion**

The purpose of the bridge pier scour design formulae are to determine the required depth of the foundation before the bridge is built, or to estimate bridge safety margins during its design life, with the goal of identifying the need for maintenance or scour protection measures' installment. Numerous pier scour formulae are available for scour depth estimation—however, they do not provide reliable results as most of them are developed from limited datasets, experimentally acquired or measured directly in casespecific conditions. Formulae developed experimentally may result in under- or overestimation of scour depth if applicated at bridge piers in prototype scale [77,133]. Thus, formula must be cautiously selected based on conditions that are similar to those for which the formula is developed [17]. In order to find out the most applicable empirical formulae, sensitivity analysis is often required [134].

Any kind of obstruction in the river, such as a bridge pier, inevitably leads to morphodynamic changes, especially in sand-bed rivers with median particle size that can be easily moved during floods [135]. Since the occurrence of scour hole is only a matter of time, it is inevitable that some sort of scour protection will have to be installed during the design life of the bridge. Riprap scour protection using launchable stone almost exclusively requires the riprap to be placed above the riverbed level in a form of riprap sloping structure, which represents additional rigid obstacles in contact with the erodible sediment bed. On the other hand, scour countermeasures could have adverse effect if not installed properly. Therefore, the effect of the proposed scour countermeasures on flow environments should be investigated in detail to achieve design effectiveness during use [136].

The riprap sloping structure contracts the flow through the bridge opening, increasing the velocity and shear stress, resulting in additional lowering of the riverbed elevation next to piers. The installment of the riprap sloping structure has an effect on the bridge hydraulics similar to the flow constrictions resulting from other bridge-related structures abutments placed in the main channel or long embankments traversing the overbanks. Sudden contraction of the oncoming flow concentrates the flow through the bridge opening, increasing the flow velocity, shear stress and turbulence in the proximity of bridge piers. Depending on the contraction ratio and associated afflux, increased flow velocity and shear stress can significantly surpass the threshold value. Thus, the riverbed deepening will be continuous until the equilibrium is reached, undermining the bridge piers. The effect of the flow contraction is local, i.e., the extent of the morphodynamic changes is rarely evident

on the wider river reach, contributing to late identification of the potential hazard to the bridge. While contraction scour is not always present in the bridge locations, depending on the flow approach section layout, it will always accompany the riprap sloping structure, as the volume of the mounded stone presents significant flow constriction. Therefore, the effect of the constriction scour must not be neglected as an influencing factor in studies focusing on scour at riprap sloping structures [137]. For many bridges protected by riprap sloping structures, scour has been found to occur further downstream, often causing a deeper scour hole [121]. Yet, the final depth and relative position of the scour hole next to the riprap sloping structure has not been investigated. Scarce data about scour monitoring next to piers with installed riprap sloping structures do not allow us to draw conclusions regarding the potentially adverse effects of such pier protection. The actual position of the scour hole could be deflected from the pier as a result of toe undermining and subsequent gradual collapsing of the riprap stones into the erosion zone, propagating the scouring process downstream of the unprotected riverbed [61].

Research focusing on the riprap primarily addresses the riprap failure mechanisms [25,64,123,138], while formulae which would estimate the scour depth downstream of the riprap protection are lacking. Similarly, the formulae used for pier scour are not applicable for piers protected by riprap sloping structures. Considering that a number of existing bridges have been retrofitted with riprap sloping structures following scour hazard occurrence, future research should be oriented towards determining the effectiveness of such scour countermeasures as well as their adverse effects on the riverbed and adjacent structures. Similarities between riprap sloping structures and river training structures exist, highlighting the need for in-detail research, as flow overtopping and additional contraction scour can significantly reduce the bridge safety under more frequent flood events driven by climate change. Although many researchers have investigated the effects of flow events on bedload transport and scour processes in laboratory studies [83,85,139–142] (for a systemic review, see also [96]), studies with field measurements are still rare [39,40], especially those with continuously measured flow waves and associated scour depths [95,143].

Based on the literature review, this study hypothesizes the following: if a large database with a sufficient range of independent hydraulic parameters, riprap geometries, and resulting scour hole depth is established, an empirical formula that provides a reliable estimation of scour depth next to bridge piers protected by riprap sloping structure can be derived. This hypothesis can be tested using field data, experimental data, numerical simulation data, as well as hybrid modeling approaches. Field data are the most reliable, but occurrence of flood events often surpasses the timeframe available for the research. Therefore, field data can be used for calibration and verification of experimental and numerical simulations on more frequent flood events. Experimental and numerical simulations provide means to simulate flood events of longer return periods, such as 1000-year floods. However, preparation of a physical model may be both time-consuming and challenging from a scaling perspective. Therefore, a physical model can be used to investigate local turbulent flow field in the vicinity of the pier. Interaction between the complex geometry and the vortex system can be explored utilizing 3D numerical simulations. The 3D CFD model overperforms the physical model in predicting turbulent velocities and sediment scour, especially if the physical model relies on a single similarity method for scaling [57]. Compared to time-consuming laboratory experiments and field investigations, numerical simulations are a relatively inexpensive and fast way to collect a large and comprehensive dataset for analysis of the scour process. Numerical models can be established for a larger flow area influencing the bridge, providing representative flow conditions at the bridge opening for the boundary conditions of the physical model. Data should be complementary in the sense that each dataset addresses the shortcomings of the others, obtaining the relevant variables across the entire range.

The state-of-the-art review and discussion on the potential research direction are the basis of the project R3PEAT (Remote Real-time Riprap Protection Erosion AssessmenT on large rivers), aiming to contribute to the field by developing empirical formulae for scour estimation next to the bridges with installed scour countermeasures in the form of the riprap sloping structure. Currently, this field is not adequately researched, while similarities with other river training structures exist. Therefore, it can be hypothesized that reliable scour equations can be developed that would combine the influencing variables used for pier scour, contraction scour and toe scour.

## **5. Conclusions**

The bridge failures data recorded worldwide indicate that hydraulic causes are the most common causes of the bridge collapses. Therefore, scour monitoring in real-time is crucial for efficient bridge management. Maximum scour depth for the design flood must be estimated, which can be compared to scour development under specific flood events and associated risk calculated. Once the scour risk is determined in the life cycle of bridges over large rivers, piers are usually retrofitted with riprap sloping structures as scour protection measures. In the literature there are numerous empirical equations developed for different pier shapes and sizes, but the ones taking into account complex flow conditions in the bridge opening with installed riprap sloping structure as scour countermeasure are lacking. This paper provides a comprehensive and relevant review of bridge scour estimation methods for piers with riprap sloping structure installed as scour countermeasure. From the state-of-the-art review, hydrodynamic conditions characteristic for such structures are singled out, separating them from single pier equations. These are shallow and on occasion supercritical flow overtopping the structure and hydraulic jump forming on the downstream slope, undermining the toe. The contributions of this work to the research field are the following:


The presented research framework is the basis of the project R3PEAT (Remote Realtime Riprap Protection Erosion AssessmenT on large rivers), aiming to contribute to the field by developing empirical formulae for scour estimation next to the bridges with installed scour countermeasures in the form of the riprap sloping structure.

**Author Contributions:** Conceptualization, A.H. and G.G.; resources, A.H., G.G., K.P. and M.L.; writing—original draft preparation, A.H., G.G., K.P. and M.L.; writing—review and editing, A.H., G.G., K.P. and M.L.; visualization, A.H.; supervision, G.G.; project administration, G.G.; funding acquisition, G.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been funded in part by the Croatian Science Foundation under the project R3PEAT (UIP-2019-04-4046).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work has been supported in part by Croatian Science Foundation under "Young Researchers' Career Development Project—Training New Doctoral Students" (DOK-01-2020).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Glossary**

*ds* scour depth


## **References**


**Antonino D'Ippolito \* , Francesco Calomino, Giancarlo Alfonsi and Agostino Lauria**

Dipartimento di Ingegneria Civile, Università della Calabria, 87036 Arcavacata, Italy; francesco.calomino@unical.it (F.C.); giancarlo.alfonsi@unical.it (G.A.); agostino.lauria@unical.it (A.L.) **\*** Correspondence: antonino.dippolito@unical.it; Tel.: +39-0984-496550

**Abstract:** Vegetation on the banks and flooding areas of watercourses significantly affects energy losses. To take the latter into account, computational models make use of resistance coefficients based on the evaluation of bed and walls roughness besides the resistance to flow offered by vegetation. This paper, after summarizing the classical approaches based on descriptions and pictures, considers the recent advancements related to the analytical methods relative both to rigid and flexible vegetation. In particular, emergent rigid vegetation is first analyzed by focusing on the methods for determining the drag coefficient, then submerged rigid vegetation is analyzed, highlighting briefly the principles on which the different models are based and recalling the comparisons made in the literature. Then, the models used in the case of both emergent and submerged rigid vegetation are highlighted. As to flexible vegetation, the paper reminds first the flow conditions that cause the vegetation to lay on the channel bed, and then the classical resistance laws that were developed for the design of irrigation canals. The most recent developments in the case of submerged and emergent flexible vegetation are then presented. Since turbulence studies should be considered as the basis of flow resistance, even though the path toward practical use is still long, the new developments in the field of 3D numerical methods are briefly reviewed, presently used to assess the characteristics of turbulence and the transport of sediments and pollutants. The use of remote sensing to map riparian vegetation and estimating biomechanical parameters is briefly analyzed. Finally, some applications are presented, aimed at highlighting, in real cases, the influence exerted by vegetation on water depth and maintenance interventions.

**Keywords:** river hydraulics; vegetation; flow resistance; turbulence; numerical methods

## **1. Introduction**

As it is well known, vegetation is an important issue on the viewpoint of catchment hydrology [1–3], since rain drops interception, evapotranspiration, and infiltration are elements to consider in surface and sub-surface water balance.

Moreover, riparian vegetation plays a key role both on the ecologic and habitat viewpoints, as well as a source of biodiversity. Indeed, vegetation creates micro-environments that can host birds, small mammals, and insects helpful to agricultural purposes, prevents fertilizers and pollutants from getting to the watercourses [4] and, because of the effect on landscape, has a significant recreational function.

On a more strictly technical viewpoint, riparian vegetation interacts with water flow, with effects both on the bank stability and on the river hydraulics. Vegetation acts on bank stability since it mechanically strengthens the soil because of the presence of roots [5–8]; moreover, it reduces the soil water content because of evapotranspiration, with the consequence of reducing interstitial pressures [9].

As for river hydraulics, vegetation clutters up part of the river cross-section [7,10,11], increases the roughness, and reduces the velocity; all this results in increased water levels and reduced water conveyance. Moreover, while the smaller average velocity on one hand

**Citation:** D'Ippolito, A.; Calomino, F.; Alfonsi, G.; Lauria, A. Flow Resistance in Open Channel Due to Vegetation at Reach Scale: A Review. *Water* **2021**, *13*, 116. https://doi.org/ 10.3390/w13020116

Received: 19 November 2020 Accepted: 29 December 2020 Published: 6 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

reduces the erosion of the riverbed and banks, on the other one, it increases the sediment deposition, what makes the water cross-sections smaller and raises flooding risk. On the scale of the hydrographic network, the general velocity reduction influences the travel time of water particles, making the peak flow control easier [12,13].

Therefore, one cannot know the general effect of vegetation in advance, but every case should be considered singularly, using proper procedures. Indeed, this effect depends on both the hydraulic and mechanical properties of the water cross-section, as of the present vegetation, that may be different according to species, phenological stage, age and, possibly, maintenance.

Numerous studies are presented in the literature on the experimental, theoretical, and computational points of view [4,14–19]. The major difficulty lies in the impossibility of studying the problem under common conditions or to draw conclusions of a general value from the experiments or from the case studies. The recent progresses in the numerical solution of the flow equations make it possible to treat single problems, but it remains difficult to adapt the codes to different conditions.

Usually, in the literature, the vegetation is considered as rigid or flexible, and according to the water level, as emergent or submerged. Flexible vegetation refers to grass, reeds, and shrubs, or, when speaking about trees, to the branch and leaf system. Combinations of the above categories, really found in natural streams and channels, are still difficult to treat.

One should note that, commonly, the river cross-sections present variable roughness along the wetted perimeter, and a typical example is given by vegetated cross-sections; moreover, one river station can be formed by more sub cross-sections, differently shaped. When using one-dimensional models of water flow, one needs a roughness coefficient representative enough that it can be computed as a weighted average of local roughness coefficients. To this purpose, several equations are present in [20] on the basis of the assumptions made on how a particular variable (discharge, velocity, contour shear stress) in the subsections is related to that in the total section.

Although the subject of this review is the flow resistance due to vegetation, however, it should be noted that research on the interaction between vegetation and flow is currently focusing on a more correct assessment of the action exerted by the shear stress on the bed and the banks [21,22], on velocity distribution [23,24], on sediment transport [25–31], on finite-sized vegetation patches [32–39], on the interaction between jets and vegetation [40,41], on processes of transport and dispersion [42–44], on evolution of patches of vegetation [45], and on one-line emergent vegetation [46].

Correctly evaluating the resistance to flow is then a major aspect of not only the river safety, regulation, and maintenance, but also of flood model calibration and validation. In the following, we will present the methods found in the literature, allowing estimation of flow resistance coefficients to input into models for flood simulation, based on different types of vegetation in the river banks and floodplains. We will then analyze briefly the three-dimensional modeling of free surface flow in the presence of vegetation and the use of remote sensing for mapping riparian vegetation and estimating biomechanical parameters. Finally, the issue of vegetation management will be dealt with, presenting a number of points of view when cutting or pruning of plants is required for safety or maintenance reasons.

## **2. Flow Resistance Equations**

According to Chow [47], the resistance to flow in artificial channels and watercourses is influenced by several factors, i.e., size and shape of the grains of the material forming the wetted perimeter, vegetation, channel irregularity, channel alignment, silting and scouring, obstruction, size and shape of channel, stage and discharge, seasonal change, suspended material, and bed load. More classifications exist, among which the best known are those by Rouse [48] and Yen [20].

As is well known, the resistance to flow can be expressed by the Darcy-Weisbach *f* friction factor, the Chézy's *C*, or the Gauckler-Strickler *k* velocity coefficients and the Manning *n* roughness coefficient; the relation among these coefficients is the following:

$$\sqrt{\frac{8}{f}} = \frac{R^{1/6}}{n\sqrt{\mathfrak{g}}} = \frac{kR^{1/6}}{\sqrt{\mathfrak{g}}} = \frac{\mathbb{C}}{\sqrt{\mathfrak{g}}} = \frac{V}{\sqrt{\mathfrak{g}R}}\tag{1}$$

where *V* is the mean flow velocity, *R* the hydraulic radius, *J* the energy line slope, and *g* the gravity acceleration.

Cowan [49], considering the case of fixed bed and ignoring the effects of the suspended solid flux on the turbulence, believes that the resistance to flow in natural water-courses depends on one hand on the shape, the dimension, and arrangement of the elements that form the roughness; and on the other hand, it depends on further dissipative effects, caused by macro-vortices produced by the flow separation in the abrupt changes of direction, crosssection shape, and vegetation; and that the overall Manning coefficient can be expressed as the sum of the relevant values.

The contribution of vegetation to the roughness coefficient can be evaluated by means of descriptive approaches, photographic comparison approaches, or by analytical methods. Even though the descriptive and photographic approaches should be considered empirical in nature, nevertheless they add to the knowledge of the river conditions and in many cases are a good way to assess a value of a friction factor or roughness/velocity coefficient.

#### **3. Descriptive and Photographic Comparison Approaches**

In these methods, one roughness coefficient is chosen on the basis of the class to which the river reach belongs. Among the descriptive methods, the best known is Chow's [47]. The author gives, for every class of channels, the minimum, average, and maximum values of Manning *n* coefficient, warning that when the channel is artificial, the average values should be used in case of good maintenance only. In Tables 5 and 6 of Chow's book ([47], p. 110–113), one can observe that the Manning coefficient is 0.018 sm−1/3 in case of the excavated channel, straight, clean, uniform cross-section with no vegetation, and 0.035 sm−1/3 in case of dense weeds. In natural streams, its values are 0.030 sm−1/3 when the cross-section is clean, and 0.045 sm−1/3 in case of weeds. We should underline that for banks with many trees, the normal values of the Manning coefficient vary in a rather large field, between 0.040 sm−1/3 and 0.150 sm<sup>−</sup>1/3.

The photographic comparison approach consists of evaluating the Manning coefficient of a given river reach on the basis of similarity with the pictures of other similar cases, for which the coefficient was estimated in ordinary or flood conditions. Chow [47] shows some pictures that can be considered as a first example, representing 24 cases concerning typical channels and one natural river, with a short description of the channel conditions and the corresponding value of Manning coefficients.

Some more manuals show color pictures carrying data and description of rivers in the USA and Australia [50–53], with features very different from each other and Manning coefficients ranging in a wide field. In addition to a qualitative description of the channel, there is some information about cross-sections and hydraulic characteristics. Almost all the watercourses exhibit riparian vegetation over the banks, but estimates of the Manning coefficients very often refer to the main channel only; indeed, the values pertaining to flooding cases are limited.

Arcement and Schneider [54] show the pictures of 15 areas subject to flooding and densely vegetated, for which roughness coefficients are evaluated. This is one of the few works trying to define the contribution of vegetation to the total flow resistance. Starting from the Cowan [49] approach and taking into account the method of vegetation density proposed by Petryk and Bosmajian [55], the authors make use of density measurements and of an effective drag coefficient, to compute the total roughness and the different

contribution to it. The contribution due to vegetation ranges between percentages of the total roughness from 64% to 81%, that is very significant values.

Some researchers think that the above methods are probably more reliable than the analytical, since, in field observations, the heterogeneity is accounted for. Nevertheless, pictures of vegetated channels are more or less limited, and, consequently, choosing a reference channel is often difficult, not to say impossible.

## **4. Analytical Methods**

In the analytical methods, plants are generally described by a characteristic diameter, *D*, and height, *hv*, or by a characteristic area of plants, *Ac*. These methods are mainly suitable in laboratory experiments, where reference is made to both natural and artificial vegetation, the latter represented by cylinders of different materials or strips of plastic. Much more limited, in this approach, are the field experiments.

In laboratory experiments, the elements representing vegetation are often distributed with parallel (called linear by some authors) or staggered patterns; there are cases where the distribution is random, as vegetation usually is in nature (Figure 1).

**Figure 1.** Plan view of (**a**) parallel, (**b**) staggered, and (**c**) random patterns.

In the case of rigid vegetation represented by cylinders, the density of vegetation *λ* is often expressed as

$$
\lambda = \frac{m\pi D^2}{4} \tag{2}
$$

where *m* is the number of cylinders per unit bed area; it is also used as a density measure of the projected plant area per unit volume, *a*, with

$$a = mD \tag{3}$$

as well as the ratio of stem diameter to spacing between the stems *D/s*.

## *4.1. Rigid Vegetation*

When the vegetation is made up by trees, in an analogy of the resistance to flow due to immersed bodies, the roughness coefficient is expressed as a function of the drag force exerted by the flow on the body, depending then on the number of trees, their arrangement, the diameters of their trunks and, where appropriate, the branch system. Usually, in the laboratory tests, the rigid vegetation is simulated by cylinders, like in Figure 2, which is practically correct when the flow does not touch the leafage.

**Figure 2.** Laboratory flume with arrays of cylinders representing rigid vegetation, University of Calabria.

4.1.1. Emergent Rigid Vegetation

In case of one isolated vertical cylinder whose axis is orthogonal to the flow direction, the resistance to flow is expressed by the drag force *FD*, computed as

$$F\_D = \frac{\rho \mathbb{C}\_D h D V^2}{2} \tag{4}$$

where *ρ* is the water density, *h* is the depth of the immersed part of the cylinder, *V* is the approach velocity, and *CD* is a drag coefficient. *CD* is a function of a stem Reynolds number computed by means of the approach velocity *V* and the cylinder diameter *D*, *ReD* = *VD <sup>ν</sup>* , ν being the water kinematic viscosity.

When the cylinder is a part of a group of elements (see Figure 3), one cannot ignore that the longitudinal and transversal interference make considerably more difficult the study of the resistance to flow [56–58]. In Figure 3, *x* is the streamwise coordinate, *z* is the vertical coordinate above the river bed, *uz* is the local time-averaged velocity, *u* is the mean velocity along the vertical, *h* is the water depth, and *hv* the vegetation height.

**Figure 3.** Side view of emergent vegetation and velocity profile.

Liu et al. [59], based on an experimental survey, described in detail the flow characteristics through rigid vegetation both in emergent and submerged conditions; they considered a staggered and a linear arrangement and for the latter referred to different densities; they also considered bottoms with different roughness. Liu et al. [59] measured the horizontal component of the velocity on several verticals between two dowels and also in the free field.

With reference to the emergent vegetation in the verticals immediately downstream of the dowels, the lower velocities are obtained; the velocities are progressively increasing, proceeding downstream and approaching the next dowel. In the free field there are higher velocities. Along each vertical, the velocities are practically constant for most of the height and gradually increase as they approach the free surface. In the case of verticals placed

between two dowels, the velocity distribution near the bed has a spike which is particularly pronounced in the vertical immediately downstream of the dowel and which attenuates in the flow direction.

The velocity spike is probably caused by the horseshoe vortex that forms at the base of the dowel, attracting the liquid of the surrounding region, which is faster, towards the base of the dowel. The fluid masses near the bottom, characterized by high speed values, mix with the higher, characterized by low speeds, creating vortices that rotate counterclockwise. These vortices, moving downstream, bring the fluid masses upwards, and this determines, as the abscissa increases, the reduction of the velocity spike. The staggered distributions determine a greater resistance than the linear [56,59].

The tests carried out with a rough bed have shown a significant change in the velocity distribution only as regards the profile immediately downstream of the dowel where there have been reductions in velocity from 30 to 130%, and also an increase of about 20% for the velocity spike near the channel bed. The measurements of the vertical component of the velocity have shown, in the vertical immediately downstream to the dowel and in the vicinity of the bed, a relatively large value which is attenuated by proceeding in the flow direction. In the vertical downstream to the dowel, Liu et al. [59] found weak downward velocity components.

As to the turbulence intensity relative to the streamwise component of the velocity, Liu et al. [59] found almost constant values along each vertical, however, its magnitude varied significantly with the location, presenting the highest values in the verticals immediately downstream to the dowels and the smaller ones in the free stream region. The turbulence generated by the dowels has a length scale much smaller than the shear generated turbulence and is quickly dissipated. In the free stream region, the turbulence intensity increases near the bed due to the shear.

To find the values of *CD* in case of sparse emergent arrays, both experimental tests and numerical simulations have been carried out, with the vegetation arrangements defined before as linear, or staggered, or random. Among the first studies, we will cite only Petryk [60], Li and Shen [56], and Petryk and Bosmajian [55].

Petryk and Bosmajian [55], to determine the Manning coefficient in a vegetated channel, implement the momentum equation for a reach, imposing that the vector sum of the weight of the control volume be equal to zero, projected on the bed direction, plus the contour resistance, plus that opposed by the tree trunk; they conclude with defining the overall Manning coefficient *n* as a function of the value relative to the soil, *nb*, and the one relative to vegetation drag coefficient *CD*, by writing

$$n = n\_b \sqrt{1 + \frac{C\_D \sum A\_i}{2gAL} \left(\frac{1}{n\_b}\right)^2 R^{4/3}}\tag{5}$$

where *L* is the reach length, *A* the area of the water cross-section, Σ*Ai* the area opposed by the vegetation to the flow. The authors consider the drag coefficient *CD* equal to 1.

Stone and Shen [61] developed a method for predicting the apparent channel velocity and the velocities in the stem layer in both emergent and submerged conditions. Four staggered arrangements of stems with varying diameter and density were used. The drag coefficient values compare well with that for a single cylinder, and the authors suggest an average value of 1.05 applicable with the average velocity in the constricted cross-section.

Baptist et al. [62], based on the balance of horizontal momentum in uniform steady flow condition, calculated the Chézy coefficient in the presence of emergent vegetation, *Ck*, by writing

$$\mathcal{C}\_{k} = \sqrt{\frac{1}{\left(1/\mathcal{C}\_{b}^{2}\right) + \left(\mathcal{C}\_{D} m Dh/2g\right)}}\tag{6}$$

where *Cb* is the Chézy coefficient of the bed. Baptist et al. [62] consider *CD = 1*.

In Table 1 are reported different expressions for the drag coefficient. In the case of Ishikawa et al. [63], Kothiary et al. [64], and D'Ippolito et al. [65] they are based on direct

measurements of the drag force. In Kothyari et al. [64], the cylinders are distributed on a grid with a staggered pattern forming angles of 30◦ with the flow direction, and the action is on one cylinder, while in the case of Ishikawa et al. [63], the angle with the flow direction was 45◦ and the action was calculated as mean on seven or thirteen cylinders, while instead in D'Ippolito et al. [65], the cylinders are in a linear arrangement and the action was calculated as mean on two to twenty five cylinders.

**Table 1.** Equations for estimating the drag coefficient *CD* in case of arrays of emergent cylinders.


(*i* is the bed slope, *α0E* and α*1E* are the Ergun coefficients, *J* is the energy slope, *Vv* is the average pore velocity *Vv = V/(1*−*λ), rv* is the vegetation-related hydraulic radius *rv* = *<sup>π</sup>* 4 (1−*λ*) *<sup>λ</sup> D*, *rv*<sup>∗</sup> is the dimensionless vegetation-related hydraulic radius *rv*<sup>∗</sup> = *g J ν*2 1/3 *rv*, *Rev* is the vegetation Reynolds number *Rev* = *Vvrv <sup>ν</sup>* , *ReD\** is the stem Reynolds number calculated with the average pore velocity *ReD\* = VvD/ν*).

> In Ishikawa et al.'s tests [63], the drag coefficients based on the experimental results differ significantly for the same stem Reynolds number, although the dependence is unclear. The authors give three equations for the drag coefficient as a function of the vegetation density, depending on the flume slope.

> Tanino and Nepf [66], on the basis of Ergun's [70] studies, proposed the equation reported in Table 1, where *α*0*<sup>E</sup>* is relative to the contribution of viscous forces on the cylinder surface, and *α*1*<sup>E</sup>* to the contribution of inertial forces deriving by the pressure drop downstream to the cylinders. The authors find that *α*1*<sup>E</sup>* is linearly varying with the density, whereas *α*0*<sup>E</sup>* is independent of the cylinder array characteristics for 0.15 ≤ *a* ≤ 0.35. These results are confirmed by Tinoco and Cowen [71] for *ReD* > 1000.

> In the tests performed by Kothyari et al. [64], the density *λ* is defined as Equation (2) and is ranging between 0.0022 and 0.0885. The mean flow velocity was estimated as the flowrate (*Q*) divided by the flume cross section and (*1* − *λ*), obtaining the so-called pore velocity *Vv*:

$$V\_v = \frac{\mathcal{Q}}{A(1-\lambda)}'\tag{7}$$

In the case of subcritical flows, the authors obtained the equation reported in Table 1 in which *CD* remarkably increases with *λ* and varies weakly with the stem Reynolds number calculated with the average pore velocity. The logarithmic increase of *CD* with *λ* ensures that, when *λ* is small, *CD* increases rapidly with *λ,* whereas, as *λ* becomes large, *CD* tends to a constant value.

Cheng and Nguyen [67], when the wall and bottom effects are negligible, introduced the vegetation-related hydraulic radius *rv* = *<sup>π</sup>* 4 (1−*λ*) *<sup>λ</sup> D*, which is a function of the vegetation density and diameter only. This vegetation-related hydraulic radius is used together with the pore velocity to define a new Reynolds number, the vegetation Reynolds number *Rev* = *Vvrv <sup>ν</sup>* . Using experimental data from several authors (random, staggered,

only two cases linear), suitably unified, they showed that the drag coefficient decreases monotonically with the increase in the vegetation Reynolds number and propose the two equations reported in Table 1, respectively, function of *Rev* and *rv\** with *rv\** dimensionless vegetation-related hydraulic radius.

Wang et al. [68], in a study for incipient bed shear stress partition in mobile bed channels, investigated the vegetation drag coefficient. By ignoring the bed surface shear stress, an empirical formula was developed by data fitting in which the drag coefficient is a function of the Reynolds number, calculated with the vegetation-related hydraulic radius of Cheng and Nguyen [67], the ratio between vegetation diameter and flow depth (*D/h*), and the vegetation density *λ*. The proposed formula is reported in Table 1.

Sonnenwald et al. [69] based on the data of Ben Meftah and Mossa [23], Stoesser et al. [72], Tanino and Nepf [66], and Tinoco and Cowen [71] proposed the equation reported in Table 1, where the coefficient of the *D* terms must have units m-1 to retain non-dimensionality.

D'Ippolito et al. [65], on the basis of 70 tests with emergent stems in a linear arrangement, proposed an equation in which *CD* is a function of the density *λ* only, valid in the field *λ* from 0.003 to 0.05, bed slope *i* from 0.48% to 2.02%, and *ReD* from 1000 to 10.000.

Figure 4 shows how the drag coefficient varies with density, and Reynolds numbers, for some of the formulas reported in Table 1. The values obtained by D'Ippolito et al. [65] are smaller with respect to those of other authors with the same *λ*, because of the different rod arrangements (square mesh against triangular mesh).

**Figure 4.** Emergent rigid vegetation—*CD* as a function of *λ* and *ReD*.

As it is easily seen, the *CD* values most frequently range between 0.5 and 2.0. One must notice that a comparison among the above equations is difficult, because in the equations of Kothyari et al. [64], *CD* has the form *CD* = *CD*(*λ*, *ReD*∗), while the equations of Cheng and Nguyen [67] have the forms *CD* = *CD*(*Rev*) and *CD* = *CD*(*rv*∗); the equation of Wang et al. [68] has the form *CD* = *CD*(*λ*, *D*/*h*, *Rev*); in the equation of Sonnenwald et al. [69], it is *CD* = *CD*(*λ*, *ReD*∗, *D*); and in the equation of D'Ippolito et al. [65], *CD* = *CD*(*λ*).

## 4.1.2. Submerged Rigid Vegetation

In the case of submerged vegetation, Baptist et al. [62] have identified in the velocity profile along a vertical four distinct areas, even though very often the profile is schematized as only two interacting zones (two-layer approach): The vegetation layer, containing the cylindrical elements representing the vegetation, and the surface layer, above them, up to the flow surface (Figure 5). In Figure 5, *us* is the mean velocity along the vertical in the surface layer, and *uv* is the mean velocity along the vertical in the vegetated layer.

**Figure 5.** Side view of submerged vegetation and velocity profile.

The flow characteristics within a set of submerged cylinders, with the same arrangement and height of the cylinders, are similar to the case of emergent vegetation [59]. Above the cylinders, the liquid flows with a higher velocity and this, at the top of the cylinders, gives rise to an inflection point. The two coflowing streams, the upper one and the one between the cylinders, give rise to a Kelvin-Helmholtz instability, which causes the liquid to rotate clockwise, causing vortices that become larger in the downstream direction, forcing the inflection point deeper into the array. In the case of sparse vegetation, the vortex affects the entire vegetated layer, whereas in the case of dense vegetation, it affects only a layer limited to the top of the dowels [7]. The longitudinal turbulence intensities reach a maximum near the top of the dowel array. It has the largest values immediately behind the dowels and decreases in the flow direction. In the free stream region, the longitudinal turbulence intensity is lowest.

Rigid submerged vegetation has been the subject of a large number of investigations [7,61,62,73–77] and comparisons [78–83]. Some researchers provided the average velocity values in the two layers, while others derived the velocity distribution and the average values [7,73,74,76]. In the vegetation layer, the streamwise velocity is usually considered constant with the flow depth [62,76], while in the surface layer various expressions were adopted for the velocity distribution [81]: The logarithmic theory [61,73,76], the Kolmogorov theory of turbulence [75], the genetic programming [62], and the representative roughness height [77,84]. Usually, to determine the constants involved in the velocity formulas, the two distributions are required to assume the same values on the separation surface between the vegetated and the surface layer. The average velocity over the entire water depth is obtained as a combination of the velocity of the vegetated layer and that of the surface layer.

Klopstra et al. [73] derive the velocity profile in the vegetated layer starting from the turbulent shear stress, described by means of a Boussinesq-type equation, and adopt the logarithmic law for the surface layer. The different constants within the model are found on the basis of three conditions at the interface (continuity of shear stress, velocity and its vertical gradient) and a condition at the bed (negligible shear stress). The value of the Chézy coefficient, derived by Klopstra et al. [73], is a rather lengthy expression, and for it we suggest referring to the original paper ([73], Equation (9)). The model has only one unknown parameter, the characteristic turbulence length scale, *αKL*, for which the following expression is proposed:

$$
\alpha\_{KL} = 0.0793 \, h\_v \ln \frac{h}{h\_v} - 0.0090 \text{ for } \alpha\_{KL} \ge 0.001 \tag{8}
$$

For this parameter, van Velzen et al. [85] instead proposed the following relationship:

$$
\alpha\_{KL} = 0.0227 \,\text{h}\_v^{0.7} \tag{9}
$$

The equations for the determination of the Chézy coefficient in the case of submerged rigid vegetation, *Cr*, proposed by Baptist et al. [62] and Cheng [84] are shown in Table 2. This table also shows the Chézy coefficient obtained from the value of the average velocity proposed by Huthoff et al. [75] and the Manning coefficient proposed by Yang and Choi [76].

**Table 2.** Equations for estimating *Cr* and *n* in case of submerged rigid vegetation.


For submerged rigid vegetation, Baptist et al. [62] derived an equation by dimensionally aware genetic programming, and this equation has been improved manually.

Yang and Choi [76] proposed a two-layer model. The momentum balance was applied to each layer and an expression for the mean velocities was proposed. The velocity was assumed uniform in the vegetation layer and logarithmic in the upper layer. In the equation reported in Table 2, *Cu* = 1 for *<sup>a</sup>* ≤ <sup>5</sup> *<sup>m</sup>*−<sup>1</sup> and *Cu* = 2 for *<sup>a</sup>* > <sup>5</sup> *<sup>m</sup>*<sup>−</sup>1.

A representative roughness height characterized by its proportionality to both the stem diameter and vegetation concentration was proposed by Cheng [84] to estimate, with the same approach of Yang and Choi [76], the average flow velocity, and thus the resistance coefficients in vegetated channels. Cheng's [84] approach was developed for submerged rigid vegetation, but also gave good results in the case of flexible vegetation.

López and García [86] employed two numerical algorithms, the *k* − *ε* and *k* − *ω* type, based on two closure equations, to model the mean flow and the turbulence structure in open-channel flows with submerged vegetation. The Manning coefficient, computed on the basis of experimental observations, shows an almost constant value close to the one corresponding to non-vegetated channels up to some threshold plant density and, once this limit is exceeded, a linear increase.

Defina and Bixio [74] extended the model by Klopstra et al. [73] and that by López and García [86] to consider both the plant geometry and drag coefficient variable with depth. A comparison between experimental data of real and artificial vegetation and the results of numerical simulations demonstrates that both models are able to reproduce the vertical profiles of velocity and shear stress within and above vegetation, whereas the turbulence characteristics are poorly predicted.

Li et al. [77] proposed the concept of the auxiliary bed and produced a dynamic twolayer model consisting of a basal layer and a suspension layer (Figure 6). The roughness of the suspension layer is defined as the product of the concentration of vegetation, *λ*, and the depth of penetration of the suspension layer in the vegetation, *δe*. The authors obtained the following expression for the Manning coefficient:

$$n = \frac{h^{5/3}}{\sqrt{\mathcal{S}}} \left[ \frac{1.96(h - h\_{\upsilon} + \delta\_{\varepsilon})^{5/3}}{\left(\lambda \delta\_{\varepsilon}\right)^{1/6}} + (1 - \lambda)(h\_{\upsilon} - \delta\_{\varepsilon})h\_{\ast}^{1/2} \right]^{-1} \tag{10}$$

where *h*<sup>∗</sup> = 2(1 − *λ*)/*CDa*, with *CD* calculated with the expression of Cheng and Nguyen [67]. For *δe*, refer to Equation (7) of Li et al. [77].

**Figure 6.** Side view of submerged vegetation with auxiliary bed and velocity profile.

Cheng [84], with reference to the flow rates, in the case of rigid submerged vegetation, compared his model with that of Stone and Shen [61], Baptist et al. [62], Huthoff et al. [75], and Yang and Choi [76]. It emerged that his model gave the best overall results while the others exhibited different behavior with varying flow rates. In particular, the Stone and Shen [61] model underestimated the flow rates especially for the higher values, while the model by Baptist et al. [62] works well for high flow rates, but yields overprediction for low flow rates, whereas Huthoff et al.'s [75] model gives the best prediction for high flow rates, but underestimates low flow rates.

Based on the Vargas-Luna et al. [79] analysis, with reference to the comparison between the measured Chézy coefficient and the one estimated with different models [61,62,73,75,76,84,85], it appears that the Klopstra et al. [73], van Velzen et al. [85], and Cheng [84] models show the best results. Although many models are based on a representation of vegetation with rigid cylinders, according to the authors, they also provide good results in the case of flexible vegetation when using the inflected vegetation height.

Pasquino and Gualtieri [81] compared the predicted and measured velocity for the Stone and Shen [61], Huthoff et al. [75], Yang and Choi [76], Cheng [84], Baptist et al. [62], and Li et al. [77] models as the density (sparse vegetation, transitional density, dense vegetation) and the degree of submergence (low submergence and high submergence) varied. They analyzed the performances of the different models on the basis of six different statistical parameters. In the case of transitional density with low submergence, all models give good results; in the case of dense vegetation with low submergence, the Huthoff et al. [75], Cheng [84], Baptist et al. [62], and Li et al. [77] models give very good results. In other cases, the number of data for estimating the statistical parameters, even if this analysis has been carried out, is less than 10.

Morri et al. [80], with the comparison between the measured and simulated mean velocities by different analytical models, highlight how the model by Huthoff et al. [75] provides the best results. It can be used for velocity lower than 0.8 m/s and underestimates the velocity values in the case of sparse vegetation, which could be explained because the bed roughness effect is neglected.

To conclude the analysis of the models proposed for submerged rigid vegetation, it must be said that they refer essentially to laboratory data and give overall good results, even if the models by Cheng [84] and Huthoff et al. [75] seem to provide the best results. One must note that comparison among the equations presented by various authors is more difficult than in the case of rigid emergent vegetation, because of the many parameters or variables involved.

#### 4.1.3. Submerged and Emergent Rigid Vegetation

Katul et al. [87] and Luhar and Nepf [11] proposed expressions to calculate the Manning coefficient for both submerged and emergent vegetation, given in Table 3.



Katul et al. [87,88] starting from the characteristics of the velocity in turbulent flow within and above rigid vegetation canopies proposed to calculate the Manning coefficient according to the relationship reported in Table 3, where *Cu* is the similarity constant (empirically estimated as 4.5) and *αKA* is the characteristic eddy size coefficient (estimated as 1 for gravel bed streams and 0.5 for canopy). The previous relationship can also be used in the case of shallow streams; more generally, it is valid for 0.2 < *h*/*D* < 7.

Luhar and Nepf [11], in agreement with Green [10], consider that the flow resistance due to aquatic vegetation depends on the blockage factor, *Bx*, which is the fraction of the channel cross-section blocked by vegetation. Additionally, for the same blockage factor *Bx*, the distribution of vegetation, in terms of distinct patches, affects the hydraulic resistance, since the interfacial area between the patches and the unobstructed flow increases when the patches are divided into smaller units. However, on the basis of observations made in natural rivers, the authors estimate velocities for the case where the blockage factor is known, but the exact distribution pattern is unknown, and it introduces up to 20% uncertainty. Luhar and Nepf [11] propose the relationships reported in Table 3 for *Bx* < 0.8, for *Bx* = 1.0 and for channels where the vegetation, of height *hv < h*, fills the entire width. In the reported equations, *kLN* = 1 *m*1/3/*s* is a constant necessary to make the equation dimensionally correct, *a* is the frontal area per unit volume, *Cf* is a non-dimensional coefficient proposed by Luhar-Nepf [11] given by *Cf* = 2*g*/*C*<sup>2</sup> where *C* is the Chezy coefficient.

## *4.2. Flexible Vegetation*

## 4.2.1. Potentially Changing Vegetation Condition

Sometimes during a flood event, as the discharge increases, the vegetation can lay over or be removed [89–92], which leads to a reduction in roughness and to an increase in the flow capacity through the section; therefore, the peak flow, which could occur later, takes lower water-surface elevations than it would have had in the case of upright vegetation. To determine under what conditions the vegetation flattens, Phillips et al. [89] referred to the stream power, defined as *SP* = *gRSV*, and to the resistance of the vegetation characterized by an index defined as the susceptibility index of the vegetation. This index is given by the product of the vegetation flexibility factor, the vegetation blocking coefficient, the vegetation distribution coefficient and, finally, the flow depth coefficient. The vegetation flexibility factor is characteristic of the type of vegetation and is equal to the bending moment determined by computing the product of the moment arm (distance from the bed to the point of application) and the force required to bend the vegetation to 45 degrees from the vertical. Phillips et al. [89], with reference to four different types of vegetation, by means of regression techniques, proposed the equations relating bending moment to vegetation height. The vegetation blocking coefficient takes into account the combined resistance associated with the vegetation and was determined by assigning a weighted value to the estimated percentage of the cross-sectional area of flow that was blocked by vegetation

for pre-flow conditions. The vegetation distribution coefficient takes into account the fact that vegetation aligned to the flow direction determines a velocities redistribution and the vegetation is affected by a smaller action, whereas in the case of uniform distribution every element is affected by the same velocity. Finally, the action of the flow to bend the vegetation depends on the water depth in relation to vegetation height. Five different categories have been identified, and the relative coefficients have been assigned according to the ratio between hydraulic radius (approximated with mean flow depth) and vegetation height. Phillips et al. [89], based on the experimental values, identified a threshold curve from which, for an assigned value of the vegetation susceptibility index, it is possible to determine the minimum value of the stream power beyond which the vegetation is layover.

Francalanci et al. [92], with reference to a river in Italy carried out the hydraulic modelling of one flood event in order to investigate the hydraulic forcing on trees, in terms of flow velocity and water depth. A conceptual model for the stability of riparian vegetation, based on the flow drag, the plant weight, and the resistance force of the root apparatus at the soil-root interface, has been used in a predictive form to investigate the range of the input variables which can promote the plant uprooting or preserve its stability.

For rigid vegetation, the action of the flow varies with the squared velocity. In case of flexible vegetation, this is not true, since the vegetation reconfigures by reducing the area projected on a plane orthogonal to the flow direction and aligning the leaves with it. The relationship between flow velocity and drag force was expressed as *FD* ∝ *V*2+*b*, where the Vogel exponent *b* [58] is a measure of the plant reconfiguration. When is *b* = −1, the drag force varies linearly with the velocity. A linear increase of drag force with the flow velocity was observed for flexible plants by direct measurement in prototype scale by Armanini et al. [93].

The friction factor due to the vegetation, *f* v, in a reach of length *L* is

$$f\_v = \frac{8F\_D h}{\rho V^2 AL} \tag{11}$$

We will first analyze the studies in the case of submerged flexible vegetation and then those related to the non-submerged vegetation.

## 4.2.2. Submerged Flexible Vegetation

The first studies on submerged flexible vegetation (Figure 7) are related to the design of irrigation canals [94]. In Figure 7, *hvf* is the bent vegetation height. Since the resistance depends on the curvature of the vegetation, the link between the Manning coefficient, *n*, and the product of velocity, *V*, and hydraulic radius, *R*, was identified on an experimental basis. This relationship depends on the type of vegetation and is practically independent on the slope of the canal and its shape. Five experimental curves have been obtained relating the Manning coefficients, also called delay coefficients, with the product VR classified as very high, high, moderate, low, or very low, and identified with the letters A, B, C, D, and E. This link was later reviewed and generalized [95–98].

**Figure 7.** Side view of submerged flexible vegetation.

The aforementioned curves, for *n* < 0.2, were interpolated with the following equation:

$$n = \frac{1}{(2.08 + 2.3x + 6\ln(10.8\,VR))}\tag{12}$$

with *x* equal to −0.5, 2, 5, 7, 11 respectively for curves A, B, C, D, and E [97].

According to Kouwen et al. [99], the friction factor in the case of submerged flexible vegetation can be represented through a semi-logarithmic function of the relative roughness, defined as the ratio of the deflected plant height, *hv f* , to the flow depth

$$\frac{\text{V}}{\text{\(\mu\_{\ast}\)}} = \text{C}\_{0} + \text{C}\_{1} \log \frac{h}{h\_{vf}} \tag{13}$$

where the coefficients *C*<sup>0</sup> and *C*<sup>1</sup> depend on *u*∗/*u*∗*<sup>c</sup>* with *u*∗*<sup>c</sup>*

$$
\mu\_{\ast c} = \min\left(0.23 \, MEI^{0.106}, \, 0.028 + 6.33MEI^2\right) \tag{14}
$$

where *M* is the stem count per unit area, *E* is the modulus of elasticity, and *I* the second moment of the cross-sectional area of the stems. The product *MEI* is defined as flexural rigidity. The deflected height depends on the drag exerted by the flow and the flexural rigidity of the vegetation. Kouwen and Unny [90] deduced the following relationship:

$$\frac{h\_{\upsilon f}}{h\_{\upsilon}} = 0.14 \left[ \frac{\left(\frac{MEI}{\gamma h \text{S}}\right)^{1/4}}{h\_{\upsilon}} \right]^{1.59} \tag{15}$$

where *γ* is the water specific weight.

Kouwen's [100] model implicitly considers the fact that the flow can cause the vegetation to lay over on the bottom of the channel without using the stream power illustrated in point 4.2.1.

Carollo et al. [91], based on the dimensional analysis, wrote the equation

$$\frac{V}{u\_\*} = \sqrt{\frac{8}{f\_v}} = \left. f\_c(m) \left(\frac{h}{h\_{vf}}\right)^{\beta\_1} \left(\frac{u\_\* h\_{vf}}{\upsilon}\right)^{\beta\_2} \left(\frac{h\_v}{h\_{vf}}\right)^{\beta\_3} \tag{16}$$

where *fc* is a function of the concentration *m* (number of stems per unit area), *hv* the vegetation height in the absence of flow, *hv f* the height of bent vegetation, and *u*<sup>∗</sup> the shear velocity. The authors carried out experimental test to determine the exponents *β*1, *β*2, *β*<sup>3</sup> and the function *fc*, which allow computation of the velocity *V*, since the bent vegetation height can be evaluated on the basis of an empirical equation.

Stephan et al. [101] analyzed, through laboratory experiments, the influence of roughness caused by aquatic vegetation on the overall flow field. The analysis of the velocity measurement shows that equivalent sand roughness and zero plane displacement of the logarithmic law are of the same order of magnitude as the mean deflected plant height.

## 4.2.3. Non-Submerged Flexible Vegetation

Studies on non-submerged flexible vegetation have been carried out, among others, by Kouwen and Fathi Moghadam [102], Freeman et al. [103], Västilä et al. [104], and Jalonen and Järvelä [105]. In particular, Kouwen and Fathi Moghadam [102], Västilä et al. [104], and Jalonen and Järvelä [105] provide an expression for the friction factor, while Freeman et al. [103] provide an expression for the Manning coefficient, given in Table 4.


**Table 4.** Equations for estimating fv and *n* in case of non-submerged flexible vegetation.

Kouwen and Fathi Moghadam [102] performed laboratory experiments on four different types of conifers both in water and in air; the authors, on the basis of dimensional analysis and a series of simplifying hypotheses, proposed to calculate the friction factor as reported in Table 4, where *ξ* takes into account the deformation of the plant, and *E* is the modulus of elasticity. The term *ξE* is called vegetation index, it is unique for each species and is obtained from the resonant frequency, mass, and length of a tree [102,106]. The authors then propose to suitably modify the resistance to flow when the projection of a plant area overlaps with that of adjacent plants.

Recently, Västilä et al. [104] proposed to calculate the friction factor, *fv*, as shown in Table 4, where *Ac* represents a characteristic area of plants, *Ab* is the bed area related to a plant, *CD<sup>χ</sup>* is a species-specific drag coefficient, *χ* depends on the area *Ac*, and *V<sup>χ</sup>* is the lowest velocity used in determining *χ*. As to the area *Ac*, Västilä et al. [104] compared three areas: The first one, *AL*, is obtained by means of the leaf area index (LAI), (defined as the ratio between the one-sided leaf area and the ground area), the second one as the area projected onto a plan orthogonal to the flow direction when the plant is under the flow action (*Ap*), and the third one as the area projected onto a plan orthogonal to the flow direction when the plant is on free air (*A*0). The results obtained by Västilä et al. [104] show that the three above characteristic areas give good results, even though the easiest area to quantify is the one relative to the leaf area index, that can be measured by laser-scanner techniques. Västilä et al. [104] determined the parameter values for Black Poplars (*CD<sup>χ</sup>* = 0.33, *V<sup>χ</sup>* = 0.1, *χ* = −1.03, 0.4 < *AL*/*Ab* < 3.21) and other foliated deciduous species of Poplars (0.43 ≤ *CD<sup>χ</sup>* ≤ 0.53, 0.10 ≤ *V<sup>χ</sup>* ≤ 0.11, −0.57 ≤ *χ* ≤ −0.9, 0.4 ≤ *AL*/*Ab* ≤ 3.2).

Jalonen and Järvelä [105], taking into account the reconfiguration properties of branches and leaves, split their contribution in the friction factor proposing the equation reported in Table 4, where the subscript *F* and *S* refer to leaves and stem, respectively. The authors give the values of the parameters (*CD<sup>χ</sup>*,*F*, *χF*, *CD<sup>χ</sup>*,*S*, *χS*) for four different species, in the velocity field from 0.1 m/s to 0.6 m/s (low) and for velocities higher than 0.6 m/s (high).

Freeman et al. [103] on the basis of dimensional analysis and extensive experimental investigation, wrote the Manning coefficient as reported in Table 4, where *A*∗ *<sup>i</sup>* is the net submerged frontal area, *As* is the total cross-sectional area of all of the stems of an individual plant, measured at a quarter of its non-deflected height, m is the density. The plant characteristic *A*∗ *<sup>i</sup>* and *As* are density-weighted average values that were measured without any bending or distortion attributable to flow. The modulus of elasticity is calculated starting from the force necessary to bend the plants at an angle of 45◦ or, more simply, by an empirical relation expressed as a function of the ratio between the height of the plant and its diameter measured at a quarter of the height.

Once again, for submerged or non-submerged flexible vegetation, comparison among the equations proposed by different authors is difficult if not impossible.

## **5. Numerical Methods**

Apart the empirical approach, turbulence study is the basis for flow resistance assessment, at least from Prandtl and von Kàrmàn theory dating back to the beginning of the last

century. Moreover, a detailed analysis of the flow fields and turbulence characteristics may be of importance as related to solid and pollutant transport.

Experimental studies of velocity profiles, like those shown in Figures 3–5, are of much interest; in some cases also the shear stress profiles are studied by means of velocity fluctuant components and Reynolds stress theory [19,107–111].

For a detailed analysis of the flow fields and turbulence characteristics, one can refer to measurements with a particle image velocimetry system (PIV) or acoustic Doppler velocimeter (ADV) probe or, furthermore, on numerical simulations, more or less detailed depending on the flow cases at hand [65,112–115]. An excellent review about the numerical models utilized for the analysis of the interaction between flow and vegetation is due to Stoesser et al. [116], and one can refer to this work for a systematic view of the different aspects of this matter. In what follows, attention is given to the more recent works on this subject.

The interaction between flow and vegetation has been studied by means of RANS (Reynolds Averaged Navier-Stokes equations [112]) and LES (Large Eddy Simulation [117]) techniques. The DNS (Direct Numerical Simulation of turbulence [118]) approach was used to analyze the interaction between fluid and cylinders, in which all the fluctuating components are computed and no closure models are used, then the swirling-strength criterion for flow-structure extraction is applied to the velocity field [119]; this method is scarcely used, mainly due to the remarkable computing resources required.

RANS models are operated on coarse grids, and the drag due to vegetation is accounted for through additional source terms in the governing equations [86,120–123]. These models require an a-priori knowledge of the drag coefficient, contain a number of empirical constants, and are fairly accurate as related to the evaluation of the mean flow field, but do not simulate the flow past single trunks or branches. Kim and Stoesser [124] demonstrated the importance of the a-priori knowledge of the drag coefficient for the correct evaluation of the resistance to flow. Their RANS simulations, carried out by adopting actual *CD* values, furnished values of the tangential shear stress on the bottom in a perfect agreement with the experimental results, while the adoption of a constant unitary value produced an underestimation of the resistance to flow as the density was greater. The coherent structures that form downstream to the cylinders or over the vegetated layer in the case of submerged vegetation can be obtained by means of the LES approach [72,117,124]. In the latter simulations, one can still refer to the resistance to flow of each single plant through the drag coefficient, or one can explicitly represent the plants by inserting in the mesh some blocked cells representing the vegetation. Stoesser et al. [72] executed a LES with reference to a companion experiment of Liu et al. [59] related to rigid emergent vegetation, obtaining interesting results. They calculated the flow characteristics with reference to reach 0.127 m long and 0.0635 m wide, with 11,604,640 computational cells. The time-averaged velocity profiles along six vertical lines resulted in a satisfactory agreement with the experimental values, both as to the streamwise and the normal components. The velocity field showed a high-velocity zone between the cylinders, the flow separation at 95◦, and a recirculation zone downstream to the cylinders with two counter-rotating vortices. Additionally, the turbulence intensities in the streamwise and vertical directions showed a good agreement with the experimental results. Then the authors executed additional simulations with reference to densities and Reynolds numbers different from the experimental case. Starting from mean velocities and pressure fields, Stoesser et al. [72] calculated the drag on the cylinders, as decomposed in pressure and shear, and the shear on the bed. The drag coefficient values resulted in a good agreement with the experimental values of Tanino and Nepf [66]. Turbulence was investigated through the instantaneous pressure fluctuations. Kim and Stoesser [124] proposed a low-resolution LES calculation in which the computing time is considerably reduced compared to a fully resolved LES. In a high-resolution LES, each vegetation element is represented by means of a curvilinear grid, while in the case of low-resolution LES, the vegetation is handled by means of a simplified immersed boundary condition on a Cartesian grid (see Figure 2 of Kim and Stoesser [124]). A square-cell grid is

preliminarily built-up. In the area interested by a generic circular obstacle, one can distinguish three types of cells, *(i)* those belonging to the obstacle for which the velocity is zero, *(ii)* those belonging to the fluid, for which the velocity is computed without any treatment, and *(iii)* those belonging to both of them (cut cells), for which the computed velocity is multiplied by the "volume fraction" (the volume of fluid divided by the total cell volume). Kim and Stoesser [124] compared the results of the low-resolution calculations with those of high-resolution calculations (different densities and Reynolds numbers of 500 and 1340 have been considered as related to a single simulation, to 534,681 and 22,994,560 grid points, respectively) and found that, in particular for low densities, the model was able to correctly represent the velocity gradients, the wakes behind the stems, and the secondary flows. In the low-resolution LES, the knowledge of the drag coefficient was not necessary, and the model could be used also to analyze complex rod configurations, namely casual distribution of the vegetation with elements very near to each other.

It should be noted that numerical simulations have been used not only for comparisons with companion experiments, but also in real situations. As an example [116,125], a RANS simulation was carried out to mirror a flood event of the Rhine river with a return period of 100 years. A reach of 3.46 km was considered, as discretized with a mesh of 258 × 64 × 12 computational cells along the streamwise, spanwise, and vertical directions, respectively, for a total of 198,144 cells. The dimension of each cell was about 13 m × 3 m × 0.5 m. The resistance to flow due to the plants was taken into account only in the momentum equation, the characteristic of the vegetation was taken into account starting from a survey related to an area of 77 m × 45 m, the Manning coefficient was calibrated on the basis of a flood event that not involved the floodplain, while the representative diameter of the vegetation was calibrated on the basis of a flood event that actually flooded the whole floodplain. The flow velocity was measured using dilution gauging techniques. The numerical results were in good agreement with the observed levels, and the values of the velocities in the zone with vegetation corresponded well to those measured by the dilution method. Additionally, the velocity distribution in the whole cross section resulted in a good agreement with the expected values. In conclusion, the turbulence study based on numerical approach is promising, even though still very far from practical use.

## **6. Hydraulic Roughness Assessment**

Riparian vegetation, in particular that located on flood areas, has very heterogeneous characteristics, both from a spatial and temporal point of view. These characteristics must be adequately included in the hydraulic–hydrological models. Conventional ground-based monitoring is often unfeasible, as these techniques are time-demanding and expensive [126], especially for large areas and when they are inaccessible. New opportunities are offered by remote sensing, which has developed considerably in recent decades and has been increasingly used in the environmental field. Some reviews have addressed the use of remote sensing in fluvial studies [127–129] and, in particular, for mapping riparian vegetation and estimating biomechanical parameters [130].

Remote sensing is based on satellite images (digital or radar) or aerial platforms (LiDAR (Light Detection and Ranging) and orto-photography). Digital satellite images, in the last decade, have reached a definition similar to those of orthophotos (the pixel sizes in the sensors Quickbird and Ikonos are equal, respectively, to 0.6 m and 1 m) making possible their application to riparian areas that, very often, are of limited size. Image classification is the process of assigning individual pixel or groups of pixels to thematic classes; it is either supervised and unsupervised.

The classification of the images is very often supervised, that is, it is based on a priori knowledge of the type of coverage, but in the case of very high resolutions, it is based on segmentation techniques [131]: Edge-finding, region-growing, knowledge-based segmentation, and a mixture of the last two.

Forzieri et al. [132] proposed a method to estimate the vegetation height and flexural rigidity for the herbaceous patterns and plant density, tree height, stem diameter, crown

base height, and crown diameter of high-forest and coppice consociations for arboreal and shrub patterns from satellite multispectral data (SPOT 5). The method is developed through four sequential steps: (1) Classification of pixel surface reflectance into five land cover classes: Mixed arboreal, shrub, herbaceous, bare soil, and water; (2) data transformation based on Principal Component Analysis of the original multispectral bands and use of only the first principal component since it explains a lot of variances; (3) identification of significant correlation structures between the main components and biomechanical properties; (4) identification/estimation/validation of the relationship (simple tri-parametric power laws) between the biomechanical properties and the normalized principal component. The vegetation hydrodynamic maps are also able to well describe the equivalent Manning's roughness coefficient as proved by comparison with simulated water stages [132].

One of the biggest limitations of optical sensors is the inability to penetrate the cloud system. Radar systems are microwave-based and do not depend on the cloud system, and are particularly useful during flooding events that usually occur in the presence of a cloud cover, allowing monitoring of the timing and spatial extent of flooding. Backscatter increases with biomass, and this makes it difficult to apply radar sensors in floodplain areas, which are usually characterized by very dense vegetation.

Satellite images provide information on the spatial variability of vegetation, but do not provide information about its vertical structure. LiDAR technology provides information on the three-dimensional structure of vegetation. Laser scanning (LS) is employed in terrestrial (TLS), airborne (ALS), and mobile (MLS) platforms. The airborne laser scanner (ALS) provides accurate information of forest canopy and ground elevations producing a digital terrain model and a digital surface model. The difference between the digital surface model and the digital terrain model gives the tree heights. Forzieri et al. [133] developed a model to identify individual tree positions, crown boundaries, and plant density using airborne LiDAR data. It needs an initial calibration phase based on a multiple attribute decision making simple additive weighting method. Jalonen et al. [134] employed multistation TLS, both in field and laboratory conditions, to derive the total plant areas of herbaceous vegetation and the vertical distribution of the total plant area of foliated woody vegetation for different levels of submergence.

## **7. Flow Resistance and Vegetation Management**

The different approaches and methods described above can be used in 1D, 2D, and 3D models.

These models can be useful in a multidisciplinary framework that includes hydrology, hydraulics, forestry, sediment transport, and ecology, which are the best vegetation management strategies. It is practically impossible to analyze the latter in the field, since one should see the flow behavior for the design discharge. To our knowledge, the only field experiments of different vegetation management is that of Errico et al. [135], who compared the effects of three different vegetation scenarios on flow velocity distribution, turbulence patterns, and flow resistance in drainage channel colonized by common reeds. The first scenario corresponded to the canopy in undisturbed conditions, the second and the third scenarios were obtained by clearing, respectively, the central part and the entire channel. The channel conveyance was obtained by clearing reeds in just the central part of the drainage channel and was comparable to that obtained by the total clearance, but with much less ecological impact and maintaining relatively high levels of turbulent intensities. Errico et al. [135] also evidenced that the most suitable models for representing natural reed canopies are those that quantify the blockage factor of the patch, rather than the effect of plant elements matrixes.

Phillips and Tadayon [90] were among the first to present some examples of vegetation management in both artificial and natural trapezoidal-shaped channels. Simulation results using HEC-RAS indicate that the design discharge for the channel for full-grown vegetation conditions would overtop channel banks and flood adjacent areas. Instead, simulations conducted for post-vegetation maintenance conditions, consisting in its partial removal,

indicate that the design discharge would remain within the channel. A common element to some of the examples shown was verifying whether the design flow has the power to lay over the bushes. If this occurred, the associated roughness component was considered negligible and not used in the determination of the Manning coefficient. They also considered a minimum amount of 30 cm of freeboard above the design water-surface elevation. The purpose of this freeboard is to mitigate risk by providing a factor of safety. As part of the various vegetation maintenance schemes or scenarios, Phillips and Tadayon [90] have highlighted how leaving the vegetation randomly distributed is aesthetically more pleasant, while leaving the trees and bushes grouped together may present a better habitat environment for wildlife.

D'Ippolito and Veltri [136] analyzed the influence exerted by vegetation on the discharge with a return time of 200 years in the mountain reach of the Crati river (Calabria, Italy). The flow profiles, determined by the use of HEC-RAS, showed that, in the reach, there is an average increase in water heights of 5% with maximum peaks of 24%. The influence exerted by herbaceous vegetation is negligible compared to that due to woody vegetation.

Luhar and Nepf [11] assert that mowing patterns that produce less interfacial area per channel length (e.g., a single continuous cut on one side of the channel) are the most effective in reducing hydraulic resistance.

Van der Sande et al. [131] created a detailed land cover map of the villages of Itteren and Borgharen in the southern part of the Netherlands by using IKONOS-2 high spatial resolution satellite imagery and combining the use of spectral, spatial, and contextual information. It is one of the first applications of land cover maps used as inputs for flood simulation models. The plant cover is classified in three classes: Natural vegetation (*n* = 0.1), deciduous forest (*n* = 0.2), and mixed forest (*n* = 0.2). Using the above maps, the model results are not very different from those obtained using less accurate maps. Van der Sande et al. [131] conclude that their use provides better results (flow direction and water depth) for less extreme events.

Abu Aly et al. [88], with reference to a gravel-cobble river from California, analyzed the effects of vegetation on velocities, depths, and extent of the flooded areas for flows ranging from 0.2 to 20 times the bankfull discharge (*QBF*). They used a two-dimensional finite volume model that solves the vertical-mediated Reynolds equations (SRH-2D) and analyzed a stretch of the water course 28.3 km long with a mesh of 1–3 m. They obtained the height of the vegetation in each cell from LiDAR data and estimated the Manning coefficient based on the approach developed by Katul et al. [87] reported above. Compared to the case of absence of vegetation, they achieved an increase in the mean water depth of 7.4% and a mean velocity decrease of 17.5% for a flow of 4 *QBF*, values that rise, respectively, to 25% and 30% for a flow rate of 22 *QBF*. The model also shows how the vegetation has a strong channelization effect on the flow, in fact the flow is diverted away from densely vegetated areas and there is an increase in the difference between mid-channel and bank velocities.

Benifei et al. [137] investigate the effects of channel and riparian vegetation on flood event, with 200-year-return-period, really happened. Since the water levels were strongly influenced by riparian vegetation on the floodplain, they analyzed the effectiveness of different flood mitigation strategies employing the model Delft3D. They use the Baptist method [62]. Based on the results of the model, Benifei et al. [137] show that the most effective flood management strategy is obtained when the high vegetation (composed by trees) inside the flooding area with a return period of 2 years is removed, avoiding the growth of the bushes (5-years-old trees) as well.

## **8. Future Research and Conclusions**

The paper is a review of a considerable number of studies on the flow resistance due to vegetation, based on experimental investigations in the presence of simulated or real vegetation; in addition, the paragraph on numerical methods also refers to phenomena on the local scale. It should be said, however, that flow resistance is a matter of detail because,

according to Tsujimoto [138], flow, sediment transport, geomorphology, and vegetation constitute an interrelated system. The analysis of this interaction requires a multidisciplinary approach that can involve a set of disciplines such as hydraulics, hydrology, sediment transport, ecology, botany, and geotechnical engineering. It should also be noted that these phenomena vary greatly over time, in fact the climate and hydrology can alter growth patterns, rates of colonization, and vegetation density. The vegetation is simultaneously a dependent variable and an independent variable. Recent research has begun to investigate the interaction between sediment and vegetation in terms of sediment transport, erosion, and deposition, as stated in the introduction, but these interactions remain poorly quantified [16] and are new broad fields of research; the same is to be said for the interaction between vegetation and river morphodynamics [139,140].

With reference to flow resistance in the case of rigid vegetation many arrangements have been investigated taking into account different variables, so that a comparison of the proposed formulas is not immediate, even though possible. The first studies started from the use of basic equations, like momentum equation, and concepts, like the drag coefficient. Later, the studies developed towards the velocity and shear stress distribution. In the case of rigid vegetation, a few variables look to be able to derive a law for flow resistance. In the case of flexible vegetation, on the other hand, the various investigations refer to specific shrub and tree varieties, and the parameters of the different formulas are related to the specific varieties investigated, as well to the measuring range of the experimental tests. The wide variety of vegetation types and hydrodynamic conditions make comparing results and drawing general conclusions that are useful in practice difficult. Shields et al. [141] present a comparison between the Manning coefficients calculated, among others, with some of the formulas presented in the previous paragraphs both in the case of rigid vegetation and flexible vegetation. From it emerges how the formulas considered *predict values within the correct order of magnitude, however, clear guidelines for selecting one approach over another are difficult to provide* [141]. These indications could be obtained, as will be better specified later, seeing the ability of the different models to reproduce real situations. For operational purposes and for the correct use of formulas, it would be useful to establish a reference catalogue in which vegetation should be classified as rigid or flexible and in which ranges of biomechanical properties (flexural rigidity): Moreover, the geometric characteristics (shape of plants, foliage density, leaf area index) of the different species present in the river bed, on the banks, and in the floodplain areas, should be reported even better identifying how these properties vary with an appropriate parameter (age, height, or other). The vegetation, both in terms of species and density, varies considerably in space and in time. The hydraulic models that can best take into account this variability are the two-dimensional models. Remote sensing can be used to define the different roughness values, which is particularly useful in large areas and/or in areas not easily reachable. To accurately estimate the main vegetation properties, the simultaneous use of different remote sensing techniques could be used. Although methodologies have been proposed to determine height, flexural rigidity, density, stem diameter, crown base height, and diameter of vegetation, this constitutes a field with broad development perspectives. Very often it is difficult to classify the different types of vegetation, and remote sensing data require calibration through direct analysis in the field.

In order to identify which of the models previously seen are the most suitable to be used in the verification and design of interventions, it would be appropriate to apply them to real cases to reproduce flood events; an example of this is that of Benifei et al. [137], even if it is limited to the use of only one model. Considering the danger of flood events, reference could be made to controlled floods or indirect measures. At the same time, the same models can be used for the field verification of different riparian vegetation management strategies [135] and, subsequently, after calibration, for the design of the interventions.

Studies should be carried out on the dynamics of flood phenomena in the sense that high discharge can cause vegetation to bend on the bottom or break it or, due to the effect of localized erosion, can determine the uprooting of plants. Consequently, peak discharge

or receding flow can have lower heights, even though this is valid only from a theoretical point of view because the material transported can accumulate on the bridge piles, causing dangerous backwater effects.

While the evaluation of flow resistance with reference to river reach is of considerable importance, it must be said that even the effect of vegetation patches exert considerable influence on flow resistance, on velocity distribution, and solid transport, and this is also a field of further research as witnessed by recent publications [34–38]. The kinetics of vegetation in terms of establishment, growth, and decay of riparian vegetation in response to dynamic hydraulic conditions are also beginning to be of interest to researchers [45,142].

While the present available results are mainly of global type, detailed numerical simulations can give insights on the characteristics of the flow fields when they interact with vegetation. This is also not a simple issue, due to the different properties that the vegetation can exhibit (rigid, submerged, etc.). Actually, the 3D numerical models can simulate very simple arrangements, and the currently required calculation resources go far beyond those used in professional practice. However, the continuous progresses in the numerical field and software development can guarantee that the detailed calculation of the flow–vegetation interaction will be soon possible. This issue involves both direct numerical simulations and numerical modeling of turbulence in this context, which may also benefit from the relevant theoretical developments and progress in observation techniques. So, in the future, they can provide more useful information for the conservation of habitats and biodiversity, maintenance of vegetation, and flooding protection.

**Author Contributions:** Conceptualization, F.C., A.D.; methodology, F.C., A.D.; writing—original draft preparation, F.C., A.D.; writing—review and editing, F.C., G.A., A.D., A.L.; supervision, F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Notation:** The following symbols are used in this paper



## **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-7662-6