*Article* **Flow Field around a Vertical Cylinder in Presence of Long Waves: An Experimental Study**

**Rosangela Basile <sup>1</sup> and Francesca De Serio 1,2,\***


**Abstract:** Long waves caused by storm surges or river floods can significantly impact marine and fluvial structures such as bridge piers. Apart from the forces that they generate on the structure, they also contribute to the formation of turbulent eddies downstream of the obstacle. This is relevant, as in this way they can affect both an erodible bottom and the ecosystem. The present study describes a medium-scale experiment, in which the propagation of two different long waves released on a steady current is investigated in the presence of a bottom-mounted rigid emergent cylinder. Velocity measurements were acquired by a Particle Image Velocimetry (PIV) system, providing instantaneous flow velocity vectors on selected 2D planes. For each experimental condition, the time-varying velocity field near the cylinder was examined in selected vertical and horizontal planes. First, we tested which analytical theory or approximated method can best represent the experimental waves. After this, we estimated the horizontal maps of velocity and vorticity downstream of the obstacle and finally processed the velocity signals by means of a wavelet-based technique, to derive the length scales of turbulent eddies. In such a way, we specifically derived how the spreading of coherent turbulent structures downstream of the cylinder depends on the features of the flume, cylinder, and wave.

**Keywords:** solitary wave; turbulent coherent structures; length scales; wavelet transform

## **1. Introduction**

In a marine and fluvial environment, civil infrastructures are often exposed to extreme wave actions. This is certainly the case of coastal and offshore facilities during storm surges, as well as of riverine bridge piers when a flood occurs. These events, typical of extreme conditions, are becoming more and more frequent and dangerous with climate change. Therefore, to provide proper structural designs and safety assessments, we must evaluate the impact of such extreme waves on these structures [1]. Generally, many experimental and numerical investigations have focused on the detection of both drag and inertia forces exerted by a flow on a vertical cylinder using the well-known Morison equation [2], requiring the detailed knowledge of the velocity field in which the cylinder is immersed.

Specifically, current–structure interaction and wave–structure interaction have been extensively examined, but separately. Among the many, Zdravkovich [3] defined the regions of disturbed flow field around a circular cylinder; Rockwell [4] characterized the spanwise flow structure around the cylinder in terms of instantaneous velocity and vorticity; Vested et al. [5] analyzed the force distribution on a vertical, circular cylinder exposed to shoaling regular waves; Antolloni et al. [6] investigated the wave-induced vortex generation around a slender, vertical cylinder; Duan et al. [7] described the shedding of vortices and the hydrodynamic forces resulting from a wave passing over a submerged circular cylinder; Chen and Wang [8] investigated with experiments and computations the wave interaction with fixed, partially submerged, vertical cylinders; Mo and Jensen [9] performed laboratory

**Citation:** Basile, R.; De Serio, F. Flow Field around a Vertical Cylinder in Presence of Long Waves: An Experimental Study. *Water* **2022**, *14*, 1945. https://doi.org/10.3390/ w14121945

Academic Editor: Giuseppe Pezzinga

Received: 27 May 2022 Accepted: 15 June 2022 Published: 17 June 2022

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

**Copyright:** © 2022 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/).

experiments and numerical tests for solitary waves breaking on a constant slope, where a vertical cylinder was installed to detect the horizontal velocity profiles. More complex structures have also been examined in presence of waves, especially to evaluate their effect in reducing the transmitted energy [10,11]. Nevertheless, limited studies have investigated the joint wave–current–structure interaction. Further, when the mutual effects of all these three components have been considered, the waves analyzed were mainly regular ones. On the contrary, based on their great lengths and periods, extreme waves are often assumed as solitary waves. Only more recent research has addressed the interaction of the solitary wave with pile structures in experimental terms and in numerical terms. As an example, Yang et al. [12] carried out a series of experiments on a solitary wave–current interaction in the presence of a vertical obstacle, focusing on the changes in wave height and water velocity. Kim et al. [13] conducted numerical simulations of a wave and current interacting with fixed offshore substructures.

It is thus clear that more investigations are needed to better understand the mechanism of the interaction between wave, current, and structure. For this reason, the present study has a twofold purpose. Firstly, we aim at contributing to this deepening of knowledge by providing a set of high-quality experimental data for solitary waves travelling on a shallow current and interfering with vertical, bottom-mounted, rigid cylindrical structures. Secondly, we aim at analyzing in detail how the flow field downstream of the obstacle is affected by the transit of the wave in terms of vorticity and turbulence. In particular, we used medium-scale laboratory models and focused on the time-varying velocity field (i) upstream of the cylinder due to the wave transit, and (ii) downstream of the cylinder due to both the wave transit and the presence of the obstruction. Specifically, the first analysis is needed to detect how the velocity distribution changes due to the contemporary presence of waves and currents, while the second is dedicated to the observation of vorticity and coherent turbulent eddies that detach and spread from the cylinder. The underlying motivation is to contribute to a thorough knowledge of the distribution of eddies downstream, characterized by different length scales. In fact, it is well-known that the condition of incipient movement of sediment particles is significantly influenced by turbulent eddies given that they can penetrate up to the channel bed [14,15]. Moreover, resembling the transport of turbulent kinetic energy, eddies can contribute to the entrainment and transport of sediments or eventually tracers, thus they can affect a mobile bottom or even the biota and the local ecosystem [16,17].

This study is arranged as follows: in the next section we introduce the experimental setup used for the investigation and the theoretical frame for the wave's description. In Section 3, first the main features of the streamwise velocities are depicted and compared with the available theories to define the best matching one based on a proper index. After this, the vertical profiles of the flow velocity are examined during the wave transit. Further, vorticity and coherent turbulent structures in the horizontal plane around the cylinder are examined. Finally in Section 4 the main conclusions are summarized.

## **2. Materials and Methods**

#### *2.1. Experimental Equipment and Procedure*

We performed the experiments at the Hydraulic Laboratory of the DICATECH of the Polytechnic University of Bari (Italy). The used rectangular flume, 25 m long and 0.4 m wide, has sidewalls and bottom fabricated of Plexiglass and is well-suited for optical measurements (Figure 1). The head tank feeding the flume can be served independently by both a low-pressure and a high-pressure water circuit. The low-pressure main circuit provided a stationary flow condition in the flume during the test (base flow). The highpressure adduction pipe released an additional water discharge in the head tank by means of an electro valve (Figure 1) managed by a process PC, equipped with a LabVIEW software. In this way, by conveniently tuning the added water release, the desired long wave was generated in the channel, superimposed to the uniform base flow. At the downstream end of the channel, a second tank is located, receiving the discharged flow. It is provided with a

triangular sharp-crested weir used to estimate the flow rate. The water level is controlled by a sloping gate at the end of the flume. In order to reduce the reflection of the generated waves, a structure with a high degree of porosity, consisting of a 2 m in length metal cage with a 0.01 m mesh filled with *d*<sup>50</sup> = 0.015 m gravel, is located in the final part of the flume on the bottom. However, the measurements of the tested waves were acquired in a time period specifically chosen to avoid any reflection.

**Figure 1.** Sketch of the channel with location of the cylinders and (*X,Z*) axes. Side view.

The physical model was designed to reproduce a possible bridge pile affected by a river flood, in Froude analogy, thus using a length scale factor equal to 1/10 (model/prototype). Therefore, two rigid cylinders having a diameter *d* = 0.02 m were located along the *Y* axis in the same transversal section, at a distance of *X* = 10.9 m from the head tank O (*X* being the longitudinal axis of symmetry of the channel as in Figure 1). They were equidistant from the *X* axis, with *Y* = 0.10 m and *Y* = −0.10 m, respectively (Figure 2).

**Figure 2.** Sketch of the longitudinal (**a**) and plan view (**b**) of the cylinders in the flume with the fixed coordinate system (*X*,*Y*,*Z*). Dashed lines define the field of views (FoV) of the PIV. Measurements not to scale.

In order to obtain the flow velocity vectors on selected 2D planes, the velocity measurements were acquired by a Particle Image Velocimetry (PIV) technique. The 2D PIV system was equipped with a FlowSense EO 4M-32 camera, a laser with pulse energy of 200 mJ at 15 Hz and a synchronizer, controlled and monitored by a PC. The system was handled in double-frame mode; the sampling frequency was settled to 8.13 Hz and the time interval between two frames of the same pair was 150 μs. Considering the importance of seeding particles to obtain good quality measurements [18], we used a neutrally buoyant silver powder as seeding tracer with a mean diameter of O(10 μm).

The data examined in the present work refer to the flow velocity (i) measured in the vertical plane (*X*,*Z*) at *Y* = 0.10 m (Figure 2a), (ii) and measured in the horizontal plane (*X*,*Y*) located at *Z* = 0.03 m from the bottom. The velocity components are: *u* along *X*, positive towards the gravel beach; *v* along *Y*, positive following *Y*; *w* along *Z*, positive upwards,

being *Z* = 0 m at the channel bottom. The field of views (FoVs) were properly selected (see dashed squares in Figure 2). The interrogation area of the images in the adaptive correlation processing was 32 × 32 pixels, thus the velocity vectors were assessed in points regularly spaced at 0.4 mm in the vertical plane and 0.8 mm in the horizontal one, resulting in a high spatial resolution. The number of such grid points for each FoV was 172 × 172, so the calibrated PIV images had the dimensions 69 mm × 69 mm and 138 mm × 138 mm, respectively, in the vertical and horizontal plane. Considering a precision in estimating the displacement within an interrogation window of 0.1 pixel, the accuracy for the detected PIV velocity was close to 1%.

The water depth in the flume was always set to be *h* = 0.1 m. The base flow rate calculated with the Thomson-type triangular weir provided a reference base flow velocity *u*<sup>0</sup> equal to 0.16 m/s for all the tests. The base flow Froude number and Reynolds number were *Fr*<sup>0</sup> = *u*0/ *gh* = 0.162 and Re0 = *u*0*B*/*υ* = 64, 000, respectively, thus representing a turbulent subcritical condition. The water kinematic viscosity *ν* was set 10−<sup>6</sup> m2/s, while *B* was the channel width.

In order to replicate a flooding wave impacting a bridge pile two different long waves were examined in the experiments named W1 and W2, respectively. Each one was generated by linearly opening and successively closing the electro-valve of the high-pressure circuit (Figure 1) for a specific temporal interval. They were characterized by a wave height *H*, wave period *T*, wavelength *L*, and an Ursell number (index of non-linearity) *Ur* = *HL*/*h*2, as shown in Table 1. The number of images acquired by the PIV was limited by technical reasons related to the storage size of the dedicated PC, thus it was set to 150 for each measurement. Consequently, the total acquisition time for each measurement was close to 18 s, which is lower than the entire wave periods. Nevertheless, it was sufficient to capture most of the ascending and descending part of each examined wave.

**Table 1.** Characteristics of the two tests.


For both W1 and W2 tests, eight identical waves were generated and successively released in the flume, each of them separated from the previous one by the time-interval necessary to reconfigure the undisturbed condition in the channel. In this way, we could operate a phase averaging technique on the measured velocity signals [17] and obtain the turbulent fluctuations, as written in the following.

Finally, it is worth noting that in our physical model the spacing between the cylinders was 10 times their diameter, and thus each cylinder behaved as isolated, not affecting the flow field around the other one. Rather, a symmetry with respect to the longitudinal *X* axis in the flow was detected. This motivated our choice to focus our current analysis only on one cylinder, specifically on cylinder 2 (Figure 2).

## *2.2. Theoretical Background for Waves Description*

To describe analytically the generated waves, based on their shape and featuring parameters, we referred to shallow waters following the classical definition [19]. In fact, in our case, the ratio between the water depth and the wavelength is *h/L* = 0.0008 for W1 and *h/L* = 0.002 for W2, i.e., less than 1/25 = 0.04.

For this scope, we used the Korteweg–de Vries equation [20] originally derived to describe shallow water waves of great wavelength and small amplitude, thus having *h* << *L*, and *H/h* << 1. Its solution provides the time-varying orbital velocity component *u* along the wave propagation direction *xi*:

$$u(\mathbf{x}\_i, t) = \frac{\omega - 4k^3(2m - 1)}{ak} + \frac{12k^2m}{a}cn^2(\mathbf{k}\mathbf{x}\_i - \omega t + \delta; m) \tag{1}$$

where *k=2π/L* is the wave number, *ω = 2π/T* is the angular frequency, *δ* is the phase and *α* is a parameter that can be scaled to any real number (but generally posed ±1 and ±6) [20]. It is called the cnoidal wave solution for it involves the Jacobi elliptic cosine function, *cn* with modulus *m* (where *m* = 0 ÷ 1). Specifically, when *m* tends to zero, Equation (1) describes the Airy linear (LIN) wave [19]. On the contrary, when *m* = 1, Equation (1) describes a solitary wave with a finite amplitude and propagating entirely above the mean water level with constant speed and without change in shape over a fairly long distance.

As an alternative solitary wave solution, we considered the regularized long-wave equation proposed by Peregrine (PER) to describe the behavior of an undular bore with a smooth wavefront and having the following expression [20]:

$$u(\mathbf{x}\_{i\prime},t) = \frac{\omega - k - 4k^2\omega}{ak} + \frac{12k\omega}{a}\text{sech}^2(k\mathbf{x}\_i - \omega t + \delta) \tag{2}$$

Moreover, the classical solitary solution of the Boussinesq equation (BOU) was applied in our study [21]:

$$u(\mathbf{x}\_{i\prime},t) = \frac{\omega^2 - c^2k^2 - 4\beta k^2}{ak^2} + \frac{12\beta k^2}{a}\text{sech}^2(k\mathbf{x}\_i - \omega t + \delta) \tag{3}$$

where typically one sets *α* = 3 and *β* = 1, while *c* is the wave celerity.

Finally, we proposed to test an engineering method not tailored on shallow waters, but widely used in practice due to its simplicity, that is the so-called linear Wheeler stretching (WH) method [22]. It attempts to correct the linear Airy theory, whose straight application is valid from sea bottom to mean water level, by stretching the vertical coordinate *Z* up to the instantaneous free surface. Thus, a transfer function between wave elevation *η* and velocity components in the time domain is provided. As for the *u* component, we have:

$$\mu(\mathbf{x}\_{i\prime}Z,t) = \frac{\pi H}{T} \frac{\cosh k[Z/(1+\eta(\mathbf{x},t)/\hbar)]}{\sinh(kh)} \cos(k\mathbf{x}\_i - \omega t) \tag{4}$$

In the present study, due to some technical limitation during the experiments, wave elevations were not assessed; therefore, we substituted in Equation (4) the exact solitary wave elevation expression computed as [21]:

$$\eta(\mathbf{x}\_i, t) = H \text{sech}^2 \left[ \frac{1}{h} \sqrt{\frac{3H}{4h}} (\mathbf{x}\_i - ct) \right] \tag{5}$$

## **3. Results and Discussion**

The results displayed in the following are based on the vector field maps assessed from pairs of PIV particle images by means of the adaptive correlation method, as implemented in data processing.

#### *3.1. Analytical Characterization of the Experimental Waves*

The models described in Section 2.2 were taken into account and their reliability in reproducing the experimental waves generated was tested. To allow a consistent comparison between theoretical and experimental waves, we specifically referred to the streamwise velocity component *u*, suitably treated. Considering a fixed position (in our case *xi* = *X*), we computed *u*(*t*) through all the aforementioned analytical models. We used the following procedure for both cases W1 and W2 to extract the velocity measurements for the wave cycle from the experimental time series. We examined the streamwise velocity measured in a large number of points distant and undisturbed by the cylinders, both in the vertical plane upstream (Figure 2a) and in the horizontal plane (Figure 2b). For each measure-

ment, the phase averaging technique was used to process the *u*(*t*) signal, so that it was decomposed in a time-averaged value *U*, a phase-averaged value *up*(*t*) and a turbulent fluctuation *u'*(*t*) [17]. The ensemble-averaged value *<u*(*t*)*>* was also taken into account with the sum < *u*(*t*) >= *U* + *up*(*t*). This phase averaging operation was applied to the eight waves measured, even if the assumption made of a solitary wave actually implies a non-periodical wave. We decided to follow this procedure because in this way (i) we could account for more robust and consistent measurements, not referring to a single episodic wave event, (ii) and we could operate statistics and extract turbulent components from the signal. Consequently, we were able to compare the trends of the theoretical *u* varying in the wave cycle with the corresponding phase-averaged experimental values *up*. As an example, for a selected point representative of the undisturbed hydrodynamic condition (far from the cylinder C2 and having coordinates *X* = 10.95 m, *Y* = 0.07 m and *Z* = 0.035 m), Figure 3 shows the comparison among the computed values of *u* and the measured *up* (normalized by the time-averaged velocity *U*). The time-averaged velocity in the two experiments is *U* = 0.28 m/s for W1 and *U* = 0.40 m/s for W2, respectively.

**Figure 3.** Comparison of computed *u* and experimental phase-averaged streamwise velocity *up* at point *X* = 10.95 m, *Y* = 0.07 m, and *Z* = 0.035 m for W1 (**a**) and W2 (**b**).

In qualitative terms, for both waves W1 and W2, the WH method seems to show the best matching with the experimental data, while the linear model provides the largest overestimations in absolute terms. To evaluate the agreement also in quantitative terms, the Wilmott index *Iw* was considered [23]:

$$Iw = 1 - \frac{\sum\_{i=1}^{N} \left(\chi\_{ci} - \chi\_i\right)^2}{\sum\_{i=1}^{N} \left(|\chi\_{ci} - \overline{\chi}| + |\chi\_i - \overline{\chi}|\right)^2} \tag{6}$$

where *χ<sup>c</sup>* and *χ* are, respectively, the computed and the experimental values; the overbar stands for the corresponding average value; and *N* is the total number of data in a wave cycle. The *Iw* index takes a value of 1 when a perfect agreement is obtained. If calculated for the whole wave cycle, it provides the highest values for the WH method for both W1 and W2 waves, resulting *Iw* = 0.99 and *Iw* = 0.97, respectively.

In any case, at a closer inspection, a worse agreement can be observed between the measured data and the WH data during the wave crest transit. Thus, we also computed *Iw* for partial data intervals, namely for the positive velocity values (typical of the wave crest transit) and for the negative velocity values (typical of the wave trough transit). The resulting indexes are shown in Table 2.


**Table 2.** Wilmott index for comparisons of the streamwise velocity trends shown in Figure 3.

It is evident that for wave W1, the WH method is the best suited one along the whole wave cycle. For wave W2, it remains the best referring to the wave trough; while for the wave crest case, the PER model better reproduces the velocities *up*. It is worth noting that we obtain analogous results (within ±2%, i.e., comparable to PIV accuracy) even when experimental velocities are examined at other undisturbed points. Finally, these results are also consistent with Iwagaki's abacus [15].

## *3.2. Upstream Vertical Profiles of the Streamwise Velocity*

For both W1 and W2 waves, we analyzed the FoVs in the vertical plane, upstream the cylinder (Figure 2a). The flow field was observed both in stationary condition (i.e., base flow) and in wave condition, during the transit of the long wave. The vertical profiles of the ensemble-averaged velocity display a quite flat and uniform trend along most *Z*, during each phase of the wave. Apart from the typical values of the ratio *h/L* and the indications of Iwagaki's abacus already written in Section 2.2, it was this observation that led us to address theoretically shallow water waves characterized by a condition of homogeneous flow along *Z*. The rapid release of the added mass into the laboratory flume head tank produced a long flooding wave that promptly affected the entire water depth from the most superficial layer towards the channel bottom. The effect of the bottom on the flow is most visible during the descending wave phase when a reduction in velocity is generally assessed over the entire depth, but is even faster in the range *Z*/*h* < 0.4. In this range, for both waves, the logarithmic trend typical of uniform flow in open channels is identifiable.

Figure 4 displays, as an example, the vertical distribution of <*u*> normalized by *U* at four different wave phases, in the profile at *X* = 10.842 m (i.e., at 0.058 m upstream of the cylinder) and *Y* = 0.10 m, for both W1 and W2. For wave W1 the phases addressed are (Figure 4a): *t*/*T* = 0.05, i.e., the wave trough; *t*/*T* = 0.3, i.e., the ascending wave phase; *t*/*T* = 0.65, i.e., the wave crest; and *t*/*T* = 0.8, i.e., the descending wave phase. Analogously, for wave W2, we have (Figure 4b): *t*/*T* = 0, *t*/*T*= 0.3, *t*/*T* = 0.55, and *t*/*T*= 0.7, respectively.

**Figure 4.** Top row: vertical profiles of the ensemble-averaged streamwise velocity <*u*>, normalized by *U*, measured at wave trough, ascending phase, crest, descending phase. Selected profile at *X* = 10.842 m and *Y* = 0.10 m. (**a**) Wave W1 and (**b**) wave W2. Bottom row (**c**): vertical profiles of <*u*> normalized by the maximum value of <*u*> detected for each profile.

The vertical profiles were obtained by stitching the instantaneous velocity maps of the two FoVs upstream of the cylinder. In fact, we originally acquired a top FoV and a bottom FoV, which overlapped 0.005 m along their common edge, as plotted in Figure 2a. In particular, for the W1 wave the lower edge of the bottom FoV was positioned at *Z* = 0.01 m, while W2 was at *Z* = 0.03 m, in order to ensure the measurement by the PIV in the upper region during wave crest transit. This transit is evident from the increasing relative heights *Z*/*h* where velocities were measured. In fact, in the trough condition the highest point investigated is close to *Z*/*h* ∼= 1 in both waves, while in the crest condition, it is close to *Z*/*h* ∼= 1.2 for W1 and *Z*/*h* ∼= 1.4 for W2 consistently with the wave heights. We observed increasing values of <*u*> due to the wave transit, as expected, reaching maxima values under the wave crest, with <*u*> = 1.24*U* for W1 and <*u*> = 1.4*U* for W2. Recalling that <*u*> is the sum of the time-averaged component and the phase-averaged velocity, we can equally write that the initial current flow at most increases by about 24% and 40%, respectively, for W1 and W2 due to the wave crest transit. In addition, approaching from the trough to the crest condition, the <*u*> values double for W1, while increasing 3.5 times for W2, meaning that a more pronounced effect of the wave on the underlying current is evident in the W2 case rather than in the W1 case.

Plotting the same vertical profiles of <*u*> normalized by the maximum value of <*u*> detected at each profile (named *max*<*u*>), we obtain Figure 4c, where all profiles tend to collapse and clearly show a flat trend for *Z*/*h* > 0.4. For *Z*/*h* < 0.4 the logarithmic trend of <*u*> is easily recognizable, especially in the W1 case, having measured points closer to

the bottom. However, for both cases a bottom boundary layer with a thickness *δ* ∼= 0.4*h* is clearly detected, regardless of the wave phase.

## *3.3. Downstream Horizontal Maps of Streamwise Velocity and Vorticity*

The horizontal maps of *<u>/U* measured at *Z* = 0.03 m are displayed in Figure 5. For the sake of brevity, only the maps corresponding to the trough and crest phases already selected and shown in Figure 4 are depicted. For greater convenience, a new local reference system (*x*,*y*) was chosen in the FoV, with the origin in the center of the cylinder.

**Figure 5.** Horizontal maps of the normalized ensemble-averaged streamwise velocity <*u*> at phases (**a**) W1 trough, (**b**) W1 crest, (**c**) W2 trough, (**d**) and W2 crest. *Z* axis is vertical and positive upward; here *Z* = 0.035 m.

Some spurious results close to the left edge of the frame are present, due to some technical limitations in the lighting system during PIV acquisition. They do not affect the signal downstream of the cylinder in the area of greatest interest, but we decided to mask them to avoid misinterpretations. For both tests, a shaded area downstream of the cylinder is always present, as expected, showing negative <*u*> values in the near wake. Due to a less intense base flow and a lower wave height in W1, the absolute values of <*u*> are lower than those in W2, both in the trough (Figure 5a) and in the crest (Figure 5b). Further, the higher *U* value for W2 determines a more extended shaded area (Figure 5c,d). In fact we observe negative values of <*u*> for *x*/*d* = 0 ÷ 3 in the W1 case and up to *x/d* ∼= 4 in the W2 case. Instead, the transversal spreading of such negative velocities is of order O(*d*) for both tests. The inversion of the sign of <*u*> delimits the shaded area, with values close to zero along its edge. As a result, a recirculation arises, with eddies detaching from the

cylinder and antisymmetric with respect to the local *X* axis, as already observed in uniform flows. However, the transit of the solitary wave, and thus the presence of a trough and a crest condition, induces a variation in the vorticity distribution and in the characteristics of coherent turbulent eddies, as explained below.

The vorticity <*Wz*> in the horizontal plane at *Z* = 0.03 m was computed starting from the ensemble-averaged streamwise <*u*> and spanwise <*v*> velocity components measured by the PIV. Figure 6 displays the corresponding horizontal maps for trough and crest phases, for W1 and W2, respectively.

**Figure 6.** Horizontal maps of the ensemble-averaged vorticity *<Wz>*, at phases (**a**) W1 trough, (**b**) W1 crest, (**c**) W2 trough, (**d**) and W2 crest. Marked points identify the locations used for the successive wavelet analysis. *Z* axis is vertical and positive upward; here *Z* = 0.035 m.

Even in these plots the areas with rough data are masked, showing unreliable results. Opposite values of vorticity are observed downstream of the cylinder, mainly clockwise (negative) for *y/d* > 0 and anticlockwise (positive) for *y/d* < 0 consistently with the stationary case of a flow investing a cylinder. During the wave cycle, from trough to crest condition, the cylinder Reynolds number Re*<sup>c</sup>* = < *u* > *d*/*υ* increases, due to the increase in <*u*>, from 3400 to 6800 for W1 and from 3600 to 11,200 for W2. Thus, a more intense fully turbulent vortex street appears, as evident by comparing Figure 6a,b as well as by comparing Figure 6c,d. The antisymmetry of vorticity along the *X* axis is generally kept, with the highest values of *Wz* located at the trailing edge of the shaded area, while a spreading of vorticity with lower intensity is noted outside. In crest conditions, the values of *Wz* increase, especially for the W2 case where they are in the range (−0.8 ÷ 0.8) Hz, while for the W1 case they remain in the range (−0.6 ÷ 0.6) Hz.

## *3.4. Turbulent Coherent Structures*

Methods based on vorticity magnitude were used fairly extensively to identify vortices in the flow. Nevertheless, it should be considered that shear layers also exhibit high vorticity magnitudes. Turbulence structures essentially evolve from the interactions between vorticity and strain rate, so in the present study, to better characterize coherent turbulent vortices in the flow, we applied the continuous wavelet transform (CWT), a power tool capable of pinpointing regions with different spatial and temporal scales in turbulent flows [24], thus providing continuous time-frequency identification of eddy structures.

Classical methods to decompose the velocity signal into frequency bands have invariably been based on the Fourier transform [25]. That is, the underlying assumption is that the same spectral components are always present in a signal. However, this approach is not particularly suited for the treatment of non-stationary signals, such as in the present case, where we distinguish and establish the length and time scales of the coherent structures residing downstream of the cylinder.

The CWT approximates a complex function as a weighted sum of simpler functions, which themselves are obtained from one simple prototype function Ψ, called the mother wavelet. Several functions can be used as the mother wavelet and we assumed the Morlet one [24,26]. In CWT, the temporal distribution of the frequency components of the signal is found by successively passing stretched and compressed versions of the function Ψ, throughout the signal. By decomposing a time series into time-frequency space, one is able to determine both the dominant modes of variability and how those modes vary in time [22].

For both waves we operated in the same way, but for convenience, we firstly focused on the results of the W2 wave. We selected six reference points in the horizontal plane (at *Z* = 0.03 m), along two different transversal sections at relative distances *x/d* = 2.5 and *x/d* = 3.5, and specifically: O1(*x/d* = 2.5, *y/d* = 1.5), O2(*x/d* = 3.5, *y/d* = 1.5), E1(*x/d* = 2.5, *y/d* = 1), E2(*x/d* = 3.5, *y/d* = 1) I1(*x/d* = 2.5, *y/d* = 0.5), and I2(*x/d* = 3.5, *y/d* = 0.5). Points O1 and O2 are outside the cylinder wake; points E1 and E2 are located at the trailing edge; points I1 and I2 are inside the wake (refer to Figure 6 for points location). For each point, the time series of the turbulent streamwise velocity *u'*(*t*) to be processed by CWT was obtained by stitching the *u'*(*t*) signals of the eight consecutive waves generated in the channel and cutting out the time interval between two succeeding waves. For the six target points, Figure 7a–f display the *u'*(*t*) time series, the corresponding scalogram obtained with CWT, and the global wavelet spectrum. The scalogram is the wavelet power spectrum where the highest values in the contour plot correspond to the most energetic frequencies of the signal. The black line represents the cone of influence for the wavelet, meaning that values external to such cone are untrusted because they are too close to the extremes of the time series. The global wavelet spectrum is obtained for each frequency scale from the integral of the wavelet power calculated over the period of investigation.

Considering points O1 and O2 (Figure 7a,d) in the wavelet power spectrum low intensities are generally observed, except than in the range of frequencies (0.0625 ÷ 0.25 s<sup>−</sup>1) where more energetic spots are evident during most of the examined time interval. They are related to the presence of large coherent structures which retain the maximum turbulent power, as can easily be seen by observing the global wavelet spectrum. According to [24,25], we can rely on the well-known Taylor's approximation to convert the time scale *1/f* into the length scale *λ*. That is *λ* = 1/ *f us*, where *us* is a proper velocity scale, greater than the turbulent velocity. In the present case, following [26], we assumed the local ensembleaveraged velocity <*u*> as *us*, resting on the assumption that it conveys downstream of the coherent structures. As a result, the above written large coherent structures outside the wake have corresponding length scales with an order of magnitude in the range O(100 ÷ <sup>10</sup>−<sup>1</sup> m), in particular varying between one-tenth of the wavelength *<sup>L</sup>*/10 and the flume width *B*, thus controlled by the flume geometry.

**Figure 7.** *Cont*.

**Figure 7.** W2 wave: signal of the turbulent streamwise velocity *u'*(*t*), spectrogram of *u'*(*t*), and global wavelet spectrum for the selected points: O1 (**a**), E1 (**b**), I1 (**c**), O2 (**d**), E2 (**e**), and I2 (**f**).

The results of the wavelet power spectrum for points E1 and E2 along the edge (Figure 7b,e) are also remarkable. Moreover, in this case, high intensities in the wavelet spectrum occur at frequencies (0.0625 ÷ 0.25 s<sup>−</sup>1), manifesting again the presence of large coherent structures similar to what was already noted for O1 and O2 points, even if characterized by lower power values. More interestingly, in E1 and E2, high power spots are visible even at higher frequencies (2 ÷ 3.5 s<sup>−</sup>1), thus, we would expect them to represent small eddies. Considering that they appear only during the peak phase of the wave (as evident by moving on the time axis in Figure 7b,e), in the computation of the associated length scales we used the <*u*> velocity assessed at the wave crest; therefore, their length scales result in the range (10−<sup>2</sup> ÷ <sup>10</sup>−<sup>1</sup> m). More precisely, the spectrograms of Figure 7b,e show that, when the wave crest is approaching, coherent structures with the size of order of the cylinder diameter O(*d*) are present along the edge; during the transit of the crest, such structures become larger and are affected by the wave height with *λ* values around O(4*H*).

Figure 7c,f display the CWT results for points I1 and I2 inside the wake but close to the edge, where a behavior analogous to what was already observed for points E1 and E2 is detected. In this case, the structures at frequencies (0.0625 ÷ 1 s−1), present during most of the recorded time period, show again a higher turbulent power, comparable to O1 (Figure 7a) and O2 (Figure 7d) cases. The power intensity of the smaller coherent structures arising at higher frequencies (2 ÷ 3.5 s<sup>−</sup>1) is quite the same as for points E1 (Figure 7b) and E2 (Figure 7e).

Referring to the coherent structures detected in the W1 case, when investigated with CWT, the turbulent power values in the spectrograms are generally halved compared with the W2 case. Nevertheless, results are similar to the W2 case in terms of frequency of occurrence of the detected coherent structures. Moreover, the order of magnitude of the length scales of the eddies is generally kept, even if their size is reduced compared with the W2 case. The reduced size of vortices behind the cylinder in the case of lower wave amplitudes was also observed by [27] for harmonic waves. For brevity, only one result for the W1 case is shown in Figure 8, referring to a representative point along the trailing edge, namely E3, at *x/d* = 3, and *y/d* = 1.

**Figure 8.** W1 wave: signal of the turbulent streamwise velocity *u'*(*t*), spectrogram of *u'*(*t*), and global wavelet spectrum for the selected point E3.

We can also consider the dimensionless Strouhal number *St* to compute the frequency of vortex shedding from the cylinder *fs*. Using the local <*u*> as the reference velocity scale and the diameter *d* as the characteristic length of the obstacle in *St* and posing *St* = 0.2, we obtain a shedding frequency varying from *fs*~2 s−<sup>1</sup> in absence of the wave (for both W2 and W1 cases) and *fs*~3 s−<sup>1</sup> for the W1 wave crest transit or *fs*~4 s−<sup>1</sup> for the W2 wave crest transit. Thus, such values are consistent with those of the eddies detected with the CWT method along the edge and inside the wake of the cylinder.

The global reading of the analysis conducted with the CWT reveals that the most energetic turbulent eddies in the horizontal plane, downstream of the cylinder, are large and have size varying in the range O(10−<sup>1</sup> ÷ <sup>10</sup><sup>0</sup> m), regardless of the position considered in the FoV. Along the edge and inside the wake other coherent structures appear at higher frequencies, consistent with the Strouhal frequency. In absence of the wave, they have the size of order of the cylinder diameter, while during the transit of the wave crest their size are of order O(4*H*), i.e., they depend on the wave height. This event, occurring during the transit of the long wave crest, is in fact due to the augmented <*u*> velocity produced by the wave itself.

## **4. Conclusions**

We experimentally investigated how the flow field can be modified, both upstream and downstream of a rigid, bottom-mounted, emergent cylinder, under the action of a travelling long wave. Two waves were tested, on a constant depth *h* = 0.1 m; that is W1 having a wave height of 0.025 m and W2 having a wave height of 0.05 m. Firstly, they were processed by means of the phase averaging technique and identified in analytical terms. In fact, referring to shallow water waves of great wavelength and small amplitude, we selected some solitary wave solutions that can describe the experimental waves, such as the Peregrine (PER) and the Boussinesq (BOU) solution of the KdV equation. We also compared both laboratory waves with the linear solution by Airy (LIN) and an approximate engineering method by Wheeler (WH). To evaluate the agreement in quantitative terms, the Wilmott index was considered. It resulted that for the W1 wave, the WH method is the best suiting one along the whole wave cycle. For the W2 wave, it provides the best matching in the wave trough description, while the PER model provides a better performance in the wave crest one.

After this characterization, the vertical profiles of the streamwise velocity were analyzed upstream of the cylinder, at selected and representative phases, i.e., wave trough, ascending wave, wave crest, and descending wave. In all cases a quite flat and homogeneous trend is observed, meaning that the long flooding waves promptly affected the entire water depth, from the most superficial layer towards the channel bottom. Regardless of the wave phase, a logarithmic trend with a thickness of 0.4*h* and identifying a bottom boundary layer is recognized. For both waves, the crest transit induces an increase in the initial current flow (by about 24% for W1 and 40% for W2).

The horizontal maps of the velocity acquired at 3 cm form the bottom were then examined, evidencing a recirculation in the near wake of the cylinder, longer for W2 than for W1, due to the higher values of the streamwise velocity. The vorticity maps in the same plane highlight the presence of opposite values downstream of the cylinder, consistently with the stationary case of a flow investing a cylinder. During the wave cycle, from trough to crest condition, a more intense fully turbulent vortex street appears. The antisymmetry of vorticity along the longitudinal axis of symmetry of the cylinder is generally kept, with the highest values of vorticity located at the trailing edge of the shaded area.

Finally, to detect possible coherent structures in the same horizontal plane, the continuous wavelet transform (CWT) was applied. It shows that large coherent structures with size varying in the range O(10−<sup>1</sup> ÷ 100 m) arise. They are controlled by the flume dimensions, regardless of the position downstream of the cylinder, and contain the highest turbulent energy. Along the edge and inside the wake, other smaller coherent structures appear with higher frequencies, consistent with the Strouhal frequency. In the absence of the wave, they

have length scales with order O(10−<sup>2</sup> m), i.e., comparable to the cylinder diameter, while during the transit of the wave crest their size increases, depending on the wave height, reaching order O(4*H*). It is therefore evident that the long wave impacting the vertical cylinder modifies the turbulence production, spreading and dissipating downstream of the cylinder. This effect must be taken into account, as it can influence connected physical phenomena, such as entrainment and transport of sediments or tracers.

**Author Contributions:** Conceptualization and methodology, F.D.S.; supervision and funding acquisition, F.D.S.; formal analysis and investigation, R.B. and F.D.S.; writing F.D.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research activity of RB was funded by the Italian Ministry of University and Research (Program PON RI Dottorati Innovativi con caratterizzazione industriale XXXIV ciclo, Decreto MIUR 2983/2018).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The technical staff of the Hydraulic Laboratory of the DICATECh is gratefully acknowledged for supporting during the research activity.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## **References**


## *Article* **Stochastic Assessment of Scour Hazard**

**David Flores-Vidriales 1, Roberto Gómez <sup>1</sup> and Dante Tolentino 2,\***


**Abstract:** Scour is the most frequent cause of bridge collapses in Mexico. Bridges located along the Mexican Pacific coast are exposed to extreme rainfall originating from tropical storms and hurricanes. Such environmental phenomena trigger sediment loss, which is known as scour. If maintenance actions are not taken after scouring events, the scour depth increases over time until the bridge collapses. A methodology to estimate the scour hazard considering both the scour–fill interaction and the Monte Carlo simulation method is proposed. The general extreme value probability distribution is used to characterize the intensity of the scouring events, the lognormal distribution is used to characterize the sedimentation process (fill), and a homogeneous Poisson process is used to forecast the occurrence of both types of events. Based on the above, several histories of scour–fill depths are made; such simulations are then used to develop time-dependent scour hazard curves. Different hazard curves associated with different time intervals are estimated for a bridge located in Oaxaca, Mexico.

**Keywords:** scour; stochastic analysis; forecasting; extreme events

## **1. Introduction**

A bridge built over a river creates an obstruction to the water flow, changing the local flow field. The erosion around the obstruction is accelerated due to the following three variables: (a) a velocity surge, (b) a difference in pressure between upstream and downstream areas produced by the perpendicular action of the flow, and (c) the generation of vortices. Scour begins when a certain level of shear stress, known as critical Shields stress, is reached. The sediment at the river bed lifts due to the action of vortices and travels downstream; such vortices produce shear stress on the sediment particles lying on the river bed and contribute to the scour process. The critical Shields stress is a function of size, shape, and material of the sediment particles.

There are two current trends that study sediment transport: trend one uses computer fluid dynamics (CFD) to solve the Navier–Stokes equations; such an approach is useful for deterministic analyses and is the best option to precisely simulate the phenomenon under study, but it demands both a great amount of computer time and a high level of capabilities. Trend two consists of a stochastic approach that is not useful to accurately forecast the outcome of a single event, but it is efficient to forecast events in probabilistic terms, and it can be a good tool to study time-dependent problems such as the case of the thickness of the sediment layer on a river bed that changes over time. The sediment transport at any given point in a time lapse, *t*, is assumed to have two possible outcomes: erosion, which is also known as scour, and accretion, which is also known as fill. Regarding the latter, a mathematical model that consists of both a data generation process (DGP) and a set of suitable simplified sediment transport equations can be used to forecast the fluctuations of the sediment layer thickness within a reasonable error.

The mechanics of sediment transport is a complex issue to study from a stochastic point of view. Thus, assumptions need to be made in order to develop a simplified mathematical approach that is useful to a probabilistic approach. Research dealing with semi-empirical

**Citation:** Flores-Vidriales, D.; Gómez, R.; Tolentino, D. Stochastic Assessment of Scour Hazard. *Water* **2022**, *14*, 273. https://doi.org/ 10.3390/w14030273

Academic Editor: Mouldi Ben Meftah

Received: 31 December 2021 Accepted: 11 January 2022 Published: 18 January 2022

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

**Copyright:** © 2022 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/).

equations for the sediment deposit is scarce. The rate at which sediment is deposited in the river bed is dependent on the availability of sediment and, thus, from the erosion of the bed and walls of the river from upstream. Variables of hydraulics, geometry, and geotechnical are fundamental for sediment transport; they also govern the behavior of secondary variables such as turbulence, flow velocity, and shear stress in the sediment particles [1]. Many studies related to the probabilistic assessment of scour hazard have been developed [2–13]. A simulation that allows forecasting the outcome of an event is a mathematical representation of events that have already happened; such simulation cannot add any new knowledge about the behavior of previous events, but it can provide useful information about possible future events.

In the last 20 years, Mexico has experienced an average of 2.5 bridge collapses per year due to scour effects [14], making scour the leading cause of bridge collapses in Mexico and most of the world [15]. The Mexican coastal highway system comprises geographic zones that are prone to extreme rainfall. For example, in 2005, hurricane Stan severely damaged highway infrastructure and caused the collapse of several bridges in the state of Chiapas. Thus, it is important to develop tools and methodologies to reduce bridge collapses due to scour.

This paper aims to propose a methodology to estimate scour hazard curves considering the following: (a) a preprocessing stage to fix missing data issues using data augmentation, (b) an information criterion for the probability distribution selection, and (c) the interaction of scour and fill events through time. The proposed methodology is illustrated in a bridge with potential scouring problems located in Mexico. Such an approach can deal with missing data, or problems such as outliers, and also has the novelty feature of addressing the accretion process that lowers maximum scour depth, making the mathematical model more complete. Moreover, the methodology is compared with the approach proposed by [5] in which the scour survival function is obtained by the SRICOS EFA method.

## **2. Data Processing**

Most databases in Mexico are somewhat flawed: they either have missing data or have outliers. Such issues introduce an additional source of uncertainty that needs to be reduced in a preprocessing stage. First, outliers are defined as any value above or below 1.5 of the interquartile range, *IQR*. Such values are removed from the data set and then treated as missing at random, as per Rubin [16]. This mechanism considers data that are not systematically related to hypothetical values of known data. Thus, it is possible to use statistical inference to obtain probable values for the missing data without having to include any information about the missing data. Hence, a methodology that produces the least amount of bias is needed, such as data augmentation, *K* nearest neighbors, multiple imputation, etc. [17–20]. In this paper, the data augmentation technique [21] is used to address the missing data issue. Such a technique has two advantages: (a) it produces a low amount of bias, and (b) it is relatively simple when compared with complex optimization techniques used for missing data analysis.

## *2.1. Preprocessing*

The available data that are recorded from the observation of a natural phenomenon such as discharge , *Q*, for a time interval *t* ∈ *Z* can be used to form an observation set, *F*, in which its members are termed events, and denoted by *ω* ∈ *F*. Let Ω be a set representing all the possible outcomes of the phenomenon. Then, the probability space of the phenomenon is (Ω, *F*, *P<sup>θ</sup>* |*θ* ∈ *I*), where the parameter space, *I*, has two sources of statistical uncertainty: one is called stochastic uncertainty, which deals with the uncertainty about a fixed parameter *θ* (Figure 1); the other is called inductive uncertainty.

**Figure 1.** Probability space, parameter space, and probability distribution.

The uncertainty associated with each event of the subset *F i*s represented by a probability function, *P*, which contains all possible outcomes of *ω* (if the sample is large enough). If the uncertainty of *P* is contained in the parameter space, then, for each *I* there is a possibility of having a different probability distribution, *Pθ*, *θ* ∈ *I*, where *θ* is a set of parameters that govern the distribution function. Hence, stochastic problems are referred to as inductive, meaning that they provide predictions about events that have not yet happened using prior observations.

## *2.2. Statistical Inference of θ*

The statistical estimation of the conditional distribution parameters of an observed data set *θ* ∈ *I* is known as inference. In order to infer the parameters for a given distribution *F* = *P<sup>θ</sup>* [*X*< *x* | *θ* ∈ *I*], an inference approach must be chosen first.

Frequentist inference is based on the idea of a limiting frequency. One of the issues of using this approach is the intrinsic epistemic uncertainty of the parameters *θ*. Such approach does not provide statements about *θ*, which is a fixed value. Thus, the uncertainty lies only in the observations *ω* (aleatory uncertainty). A frequentist approach usually has a predetermined sample size because *p* values are calculated over a sample space Ω [22]. In Bayesian inference, *θ* is considered a random variable. Some knowledge about the distribution of *θ* is assumed to be available, which is called a prior distribution. Then, after observing some data *ω*, the distribution of *θ* must be updated to match the knowledge obtained during the observation process. The Fisherian approach uses the likelihood function, *L*, to compare different values of *θ* with the probability of the observed data so that *L*(*θ*) = *P<sup>θ</sup>* [17]. Such likelihood provides a relative measure of the set of parameters *θ*, whose exact value is impossible to know. Hence, the parameter values that maximize the likelihood function *L*(*θ*) must be found. Based on the above, the Fisherian approach is used because maximizing the likelihood function *L*(*θ*) does not require previous knowledge of the phenomena, and *θ* does not have an underlying distribution to propose. The combined likelihood for a given distribution of *n* parameters can be expressed as

$$L(\theta|data) = \prod\_{n=1} P\_{\theta}(data|\theta) \tag{1}$$

The maximum of the likelihood function is obtained by two steps. The log-likelihood function is obtained, and then the score function, *S*, is calculated as follows:

$$S(\theta) \equiv \nabla \log \log \, L(\theta) \tag{2}$$

Hence, we obtain the maximum likelihood estimators, MLE, and the solution of *θ*, *S*(*θ*) = 0. MLE is not flawless, and for small data samples, the Simpson's paradox could occur. Moreover, most of the real score functions do not have a closed solution. Thus, a numerical method must be used to approach the solution. Maximum likelihood estimators

(MLE) are often the best way to deal with stochastic uncertainty because they require no prior knowledge about the phenomena.

#### *2.3. Missing Data Inference*

The second source of uncertainty is known as inductive uncertainty and deals with the lack of information of the observations available. Rubin [23] uses multiple imputation techniques to reduce the inductive uncertainty. When a small data set is used, stochastic uncertainty becomes relevant, and inductive uncertainty becomes critical [17]. Hence, it is important to rely on data sets with reasonable sizes. In order to reduce the inductive uncertainty, the missing at random mechanism (MAR) is used as follows [21]:

$$f(Q\_\prime M \mid \theta, \psi) = f(Q \mid \theta) \times f(M \mid Q\_\prime \psi \mid) \tag{3}$$

where *M* is the missing data indicator with values between 0 and 1, which represent an observed value and no observed value, respectively; *f*(*ω*, *M* |*θ*, *ψ* ) is the joint probability distribution of (*ω*, *M*). The set of parameters that governs the distribution of the observed data in the subset *Fobs* is *θ*, whereas for the missing data *Fmiss* is *ψ*. Hence, under the premise of MAR, given a set of observed data *ω* that is augmented by a quantity *Z*, the augmented posterior distribution *P*(*θ*|*ω*, *Z*) can be calculated as

$$f(M \mid F\_{\text{obs}}, F\_{\text{miss}}, \psi) = f(M \mid F\_{\text{obs}}, \psi \;) \times f(F\_{\text{obs}} \; | \theta) \tag{4}$$

which means that the parameters from the data model *θ* are independent of the parameters in the missing data. Authors in [24,25] developed the following data augmentation algorithm:

$$F\_{\rm miss}{}^{(t+1)} \approx P\_r(F\_{\rm miss}|F\_{\rm obs}, \theta^t) \tag{5}$$

$$
\theta^{(t+1)} \approx P\_r(F|Q\_{obs} \; \; \; F\_{miss}^{-(t+1)}) \tag{6}
$$

where *Pr* is the probability associated to *Fmiss*. Equation (5) is the imputation step, which generates imputed values for *Fmiss* using the missing data probability distribution, observed data *Fobs*, and the set of parameters *θ* at the given iteration *t*. Equation (6) generates parameter values for the posterior distribution, using the observed data *Fobs* and the imputed values of *Fmiss* at the iteration *t* + 1. This procedure is performed iteratively until convergence is reached, and a new "complete" data set *Ff ull* is generated.

## *2.4. Suitable Probability Distribution*

The suitable distribution function of the data *F*(ω, *t* ∈ *T*) is the backbone of the DGP. In order to select the best fitting distribution to the data set, the Akaike information criterion [26], *AIC*, is used, which selects the best probability distribution based on the parsimony principle. *AIC* is obtained using the following equation:

$$AIC = -2n\,\underline{L}(\theta|F) + 2k \tag{7}$$

where *n* is the number of observations in *F*, *L*(*θ*|*F*) is the log-likelihood function, and *k* is the number of parameters in the mathematical model.

## *2.5. Simulation of Events for the Proposed Approach*

The block maxima approach (BMA) is used to estimate the intensity and occurrence of extreme river discharges associated with extreme rain falls. Such an approach divides the observations of the "complete" data set, *Ff ull*, into non-overlapping blocks that restrict the observation of the maximum on such periods. One-year blocks are used in this research since there is no seasonality in the maximum annual discharge [27]. The generalized

extreme value distribution (GEV) is used to describe the annual maximum discharge, whose CDF is

$$F(\mathbf{x} \mid \sigma, \mu, \xi \neq 0) = \exp\left\{-\left[1 + \xi \left(\frac{\mathbf{x} - \mu}{\sigma}\right)\right]^{-\frac{1}{\xi}}\right\} \tag{8}$$

$$F(\mathbf{x} \mid \sigma, \mu, \xi = 0) = \exp\left\{-\frac{\mathbf{x} - \mu}{\sigma}\right\} \tag{9}$$

where *μ*, *σ*, *ξ* are the location, scale, and shape parameters, respectively. The law of small numbers provided the basis for the selection of the Poisson stochastic process for the time of occurrence of events that trigger scour and fill. Using this law, the stochastic process of events that exceed a given threshold is

$$P(k \ge K \mid t = T) = e^{-\lambda T} \frac{\left(\lambda T\right)^k}{k!} \tag{10}$$

where

$$
\lambda\_{\varepsilon} = \frac{\sum\_{i=1}^{n} m \mid \forall \ m \ge \psi}{n} \tag{11}
$$

and

$$
\psi \cong 0.85IQR \tag{12}
$$

Since there is only interest in the probability of occurrence of a single event [28], the cumulative distribution function (CDF) of the arrivals *Fa*(*t*) in the Homogeneous Poisson Process (HPP) is obtained as follows:

$$F\_{\mathfrak{a}}(t) = \int\_0^\infty 1 - e^{-\lambda t} \, d\_t \tag{13}$$

Having characterized the intensity of the events and its occurrence, the joint probability distribution (JPD) is

$$f\_{\mathbf{X}\mathbf{y}}(\mathbf{x},\mathbf{y}) = \sum\_{y=1}^{n} \frac{\lambda^y e^{-\lambda}}{y!} \left( \int\_{-\infty}^{\infty} \frac{1}{\sigma} \exp\left\{-\left[1 + \tilde{\xi}\left(\frac{\mathbf{x}-\mu}{\sigma}\right)\right]^{-\frac{1}{\tilde{\xi}}}\right\} \left[1 + \tilde{\xi}\_1\left(\frac{\mathbf{x}-\mu}{\sigma}\right)\right]^{-1-\frac{1}{\tilde{\xi}}} d\mathbf{x}\right) \tag{14}$$

where *y* is the number of events in a one-year block maxima.

Events that trigger fill have a random nature and have a somewhat fixed occurrence in a year. Such events occur both at the beginning and at the end of the rainy season, when the discharge is at its weakest intensity, but it is still sufficiently strong to transport sediment from the upstream river flow. Therefore, a stochastic model is needed to account for these observations that will trigger the accretion process. The intensity of the discharge, *qt*, is considered to follow a lognormal distribution function whose CDF is

$$F(\mathbf{x} \mid \mu, \sigma) = \frac{1}{\sigma \sqrt{2\pi}} \int\_0^\mathbf{x} \frac{1}{t} \exp\left\{ \frac{-\log(t-\mu)^2}{2\sigma^2} \right\} dt \tag{15}$$

where *μ* and *σ* are the shape and location parameters, respectively. The occurrence is obtained in a similar way as extreme events, with the difference in the threshold *ψ*<sup>2</sup> as follows:

$$
\lambda\_1 = \frac{\sum\_{i=1}^n m \mid \forall \ m \le \psi\_2}{n} \tag{16}
$$

where *ψ*<sup>2</sup>

$$
\psi\_2 \cong 1.2 \,\text{IQR} \tag{17}
$$

Thus, the JPD for the fill is

$$f\_{xy}(y,t) = \sum\_{y=1}^{n} \frac{\lambda^y e^{-\lambda}}{y!} (\frac{1}{\sigma \sqrt{2\pi}} \int\_0^{\infty} \frac{1}{t} \exp\left\{ \frac{-\log(t-\mu)^2}{2\sigma^2} \right\} dt \tag{18}$$

#### *2.6. Simulation of Time Series for SRICOS-EFA Method*

A time series process is a set of observations that are well defined and have regular intervals; thus, it is a set of random variables {*Xt*} indexed by integers *t* (often associated by dates). A stochastic process is a parametrized collection of random variables {*Xt*}*t*∈*<sup>T</sup>* defined on a probability space (Ω, *F*, *P<sup>θ</sup>* |*θ* ∈ *I*). Then, a time series can be understood as a realization of a stochastic process. Time series are often studied on the time domain using an autoregressive approach. According with the above, an autoregressive moving integrated average model (ARIMA) is used as follows:

$$y\_t = \varepsilon + \phi\_1 y\_{t-1} + \dots \,\_p\phi\_p y\_{t-p}' + \theta\_p \varepsilon\_{t-1} + \dots \,\_q\theta\_t \varepsilon\_{t-q} + \varepsilon\_t \tag{19}$$

where *y<sup>t</sup>* is the difference series, *p* is the order of the autoregressive, *d* is the degree of the first differencing, and *q* is the order of the moving average. The ARIMA model that best fits the data is selected by the AIC criterion [26].

Once the model is selected, simulations are performed under a random walk approach for the extreme events data set (maximum annual discharge) and then used to estimate the scour depth using the SRICOS-EFA methodology [5]. These results are then compared to the proposed approach explained in Section 2.5.

## **3. Scour and Fill**

*3.1. Scour*

Local scour depth, *ys* around a bridge pier is obtained using the one-dimensional HEC 18 approach [19] as follows:

$$y\_s = y\_n K\_4 \left[ 2K\_1 K\_2 K\_3 \left(\frac{b}{y\_n}\right)^{0.65} F^{0.43} \right] \tag{20}$$

where *yn* is the height of the stream (derived from the Manning equation), *K*<sup>4</sup> is a bed material coefficient, *K*<sup>1</sup> is a correction factor for the pier nose shape, *K*<sup>2</sup> is a correction factor for the angle of attack of the flow, *K*<sup>3</sup> is a correction factor for bed condition, *b* is the pier width, and *F* is the Froude number. The local scour area, *As*, and wetted perimeter of the scour hole, *Pmh*, are calculated as

$$A\_s \cong y\_s^{\;2} \tag{21}$$

$$P\_{mh} \approx \text{2.33} y\_s \tag{22}$$

When calculated, these values are added to the hydraulic area and wetted perimeter after each scouring event.

## *3.2. Fill*

Haschenburguer [29] proposed a model to estimate the mean fill using both the Shields stress and the hydraulic radius, which are both considered discharge-dependent variables. The Haschenburguer fill model estimates the sediment deposit from a given intensity for small-scale events in stable channels. The model uses the Shields stress as the fundamental variable, and it does not limit the amount of accretion, a fact that could lead to unrealistic heights of the sediment layer. Such issue is also present in the HEC 18 equations for scour. Hence, it is considered to have a sufficient degree of approximation [30,31].

The scour–fill process is controlled by the transfer of material in the river bed. Thus, scour depths are time-dependent random variables. The underlying process which governs sediment transport has been a constant topic for research [32–34]. A Poisson distribution

coupled with a homogeneous Poisson process (HPP) is used, so that the occurrence of the phenomenon can be expressed with a semiempirical PDF as [35]

$$f\left(y\_f\right) = \phi e^{-\phi y\_{yf}} \tag{23}$$

where *f yf* is the proportion of stream bed fill for a given flow depth, *yy f* , which is the flow depth obtained from the minimum discharge at the start of the rainy season.

$$
\phi = 3.3e^{-1.52\tau^\*/\tau} \tag{24}
$$

where the Shields mean stress, *τ*∗, is

$$\tau^\* = \frac{\rho\_{\rm w} R\_H S}{(\rho\_s - \rho\_{\rm w}) D\_{50}} \tag{25}$$

where *ρ<sup>w</sup>* is the density of the sediment, *ρ<sup>s</sup>* is the density of water at 20 ◦C, *D*<sup>50</sup> is the mean size value of the particles, and *RH* is the hydraulic ratio at the beginning and at the end of the rainy season. The hydraulic ratio is different from the one considered in the scour process, and *S* is the slope of the hydraulic grade line. A shear velocity criterion for incipient motion of sediment is used as follows [36]:

$$\tau \cong 0.215 + \frac{6.79}{D\_{50}^{1.70}} - 0.75 \times \exp\left(-2.61 \times 10^{-3} \times D\_{50}\right) \tag{26}$$

The local fill area, *As f* , and the wetted perimeter, *Pm f* , are subtracted after each event and are calculated as

$$A\_{sf} \cong f\left(y\_f\right)^2\tag{27}$$

$$P\_{mf} \approx 2.33f\left(y\_f\right) \tag{28}$$

## **4. Scour Hazard Curves**

## *4.1. Simulation of Random Events*

Both the hydraulic radius (*Ht*) and the Froude number (*F*) change through time as scour and fill develop. Hence, an analytical solution to the stochastic process is difficult to obtain. A numerical approximation to the solution is used instead. Monte Carlo simulation (MCS) is a tool for the analysis of complex systems. The MCS approach used in this research simulates a large number of independent histories of occurrence–intensity, which are then used to estimate the time progression of scour. Each simulated history is one possible outcome in the probability space. Each history behaves differently because of the intrinsic stochastic nature of the state variables (scour depth, flow, Froude number, etc.). If a sufficient number of simulated histories is known, the probability function *P* of the phenomena can be approximated (law of large numbers).

Usually scour and fill in a bridge system (a collection of variables and components) are dealt with in a discrete time fashion, especially when a probabilistic approach is used. In this research, a discrete-event system is used, assuming that the state variables change instantly through discrete points in the simulation time. Since the form of the solution is unknown, a raw sampling approach was used.

#### *4.2. Hazard Curves*

Hazard curves show the intensity of scour depth and its corresponding probability of exceedance. The probabilistic scour assessment provides the scour intensity corresponding to a target hazard level (annual probability of exceedance). The proposed hazard curve for scour is site specific, which means that there is not a variable that deals with the sourceto-site distance, which is a consequence of the fundamental approach to the phenomena. The location of the phenomenon does not have an influence on the model in the approach

proposed. Therefore, the effects of the phenomenon govern the model, and its location is negligible. The annual rate of exceedance for a given scour depth can be estimated as follows [37]:

$$P[y\_s(te) > y] = \frac{\lambda\_{\mathcal{C}}}{\tilde{\mathcal{G}}0} \left[ \tilde{\mathcal{G}}\_{yt} \right] \tag{29}$$

where *ξ*<sup>0</sup> is the mean number of events that exceed a given scour threshold, (*ψs*), during an interval *t*. *P*[*ys*(*te*) > *y*] is the probability that one or more scour depth events *ys*(*te*) are greater than a threshold (*y*), and *ξyt* is the exceedance rate given a time interval *t* and is estimated as

$$\mathcal{L}\_{yt} = 1 - (1 - P[y\_s(te) > y])^\dagger \tag{30}$$

Thus, the annual rate of exceedance can be obtained by

$$\mathfrak{F}\_0 = \frac{\left[1 - \int\_{-\infty}^{\infty} \exp\left[-\left(1 + \xi \frac{\mathbf{x} - \mu}{\sigma}\right)^{-\frac{1}{\xi}}\right] d\mathbf{x}\right] \* \mathfrak{F}\_{yt}}{\lambda\_{\varepsilon}} \tag{31}$$

#### **5. Case Study**

The bridge is located in the state of Oaxaca in one of the highest scour-prone zones in Mexico. The bridge is a four-span, simply supported bridge; each span has a length of 56 m, as shown in Figure 2a. Bridge piers and foundation piles are built with reinforced concrete in a multi-column bent type with circular columns of 1.2 m in diameter, as seen in Figure 2a. A simplified trapezoidal cross-section is used, as shown in Figure 2b.

**Figure 2.** Study case: (**a**) typical pier of the bridge; (**b**) simplified cross-section of the river.

The properties of the simplified cross-section for a given value of *H* are as follows:

$$A\_h = H\left\{b + \frac{H(m\_1 + m\_2)}{2}\right\} \tag{32}$$

$$P\_m = b + H\left[\sqrt{1 + m\_1^2} + \sqrt{1 + m\_2^2}\right] \tag{33}$$

$$Q = \frac{1}{n} A\_h R\_h^{\frac{2}{3}} S\_h^{\frac{1}{2}} \tag{34}$$

where *Ah* is the hydraulic area, *Pm* is the wetted perimeter, *Rh* = *Aht*/*Pm*, and *H* is the height of the stream for a given discharge *Q*. Equation (34) can be written as

$$Q = \frac{1}{n} \left[ \frac{\left( H \left\{ b + \frac{H(m\_1 + m\_2)}{2} \right\} \right)^{5/3}}{\left( b + H \left[ \sqrt{1 + m\_1^2} + \sqrt{1 + m\_2^2} \right] \right)^{2/3}} \right] S\_h^{\frac{1}{2}} \tag{35}$$

*Ht* is the new value of *ys*, and in a similar way *τ*∗ is deducted from the river discharge for recurring events; thus, two different hydraulic ratios are used to find the values of *ys* and *τ*∗. The values used in Equation (35) are taken from the Monte Carlo simulation. For this case study *b* = 80 *m*, *m*<sup>1</sup> = 0.8, *m*<sup>2</sup> = 1, *s* = 0.0068, and *n* = 0.025.

## *5.1. Data Augmentation*

Figure 3 shows the river discharge, *Q*, obtained from the National Superficial Water Database for this river. Both the maximum annual record and the minimum annual discharge are presented. Scour phenomenon is controlled by events that are infrequent in nature and have a high magnitude; these are known as "extreme events". Extreme events are represented in the maximum annual discharge data base, and they are used to forecast scour depths. Unlike scour, fill is controlled by recurring events at the beginning and end of the rainy season; then, the minimum discharge is used to forecast sediment deposit depths. From the data set shown in Figure 3, a "complete" data set is estimated using the data augmentation algorithm described in Section 2. Results are shown in Figure 4.

**Figure 3.** Discharge data.

**Figure 4.** Comparison between original data set and data augmentation set.

#### *5.2. Fit to a Probability Distribution*

Different probability distribution functions (PDFs) are tested to find the appropriate PDF that best characterizes the maximum annual discharge and, in a similar fashion, the minimum annual discharge. Therefore, the PDF with the lowest AIC is selected. Figure 5a shows the AIC for different PDFs such as generalized extreme value (GEV), Weibull 2P, normal, exponential, and logistic. It is noticed that the general extreme value distribution is the best fit for the maximum annual discharge. Figure 5b shows different PDFs to characterize the minimum annual discharge data at the beginning and end of the rainy season. Figure 5b indicates that the lognormal distribution presents the minimum AIC; thus, it is the best fit for the minimum flow.

The selected PDFs that characterize the maximum and minimum annual discharge are represented in probability papers. The probability paper can be defined as a mathematical– graphic technique to verify if the data follow a PDF. Figure 6a shows the maximum annual discharge and the GEV PDF with a continuous line in probability paper. Figure 6b shows the probability paper of the minimum annual discharge with a lognormal PDF.

**Figure 6.** Probability papers for (**a**) maximum annual discharge; (**b**) minimum annual discharge.

Once the probability distribution for the intensity of events has been found, the mean number of events at a specific time interval is found, as described in Section 2.5, and then the joint probability distribution (JPD) is used with the simulations, as shown in Figure 7.

**Figure 7.** Joint PDF of maximum annual flow and number of events per year.

Table 1 shows all the constants and variables needed to carry out the simulations of scour and fill events; *k*1, *k*<sup>2</sup> = 1 for the scour process.

**Table 1.** Variables and their distributions.


## *5.3. Simulation of Events*

Monte Carlo simulations were carried out to generate random events that could trigger the scour and fill process. In this context, the necessary number of simulations to have a probability *α* ≡ 0.99 providing a percentage error *ε* ≤ 5 %, is 1060, which is the number of simulations carried out to solve the stochastic processes for each simulated time interval. Figure 8a,b shows the outcome of the simulations for a time interval equal to 170 years for the extreme and recurring events. The discharge extreme events are several orders of magnitude bigger than the recurring events. The simulations of occurrences and intensity are used for the estimation of scour and fill in a discrete time model.

**Figure 8.** Simulation outcome of occurrences and intensity in accordance with both phenomena and its occurrence rates λ: (**a**) maximum annual streamflow; (**b**) minimum annual streamflow.

For comparison's sake, a set of 1024 simulations of 175 years of daily mean discharge is investigated using an ARIMA model as described in Section 2.6. The mean daily discharge data used are shown in Figure 9a. Figure 9b shows the outcome of the simulations. It can be noted that Figure 9a has outliers that have a negative impact on the forecasting precision. Since the data sets used for time series are generally much larger than that used to forecast extreme events (because of the BMA approach), outliers are more common and hard to deal with. This is a limitation on the use of time series such as ARIMA. In addition, when the data series is large and complex, advanced forecasting techniques such as wavelet multiresolution analysis should be used to improve accuracy.

## *5.4. Scour Survival Function*

Considering that the bridges in Mexico are designed for a service life of 100 years, several simulated histories were developed for time intervals of 20, 30, 40, 60, 75, 90, 110, 130, 150, and 170 years, in order to track the degradation of thickness of the sediment layer. Figure 10a shows one simulation of the thickness of the sediment layer *ys*(*t*) in the river bed as a function of time; positive values denote scour, whereas negative values denote fill. Figure 10a shows the result of a single simulation for *t* = 170 years (simulation starts at the end of the recorded data year 64). Scouring events and fill events are added in the

simulation. Since scour is controlled by extreme events, its magnitude is bigger than fill. Figure 10b shows the outcome of 1064 simulations.

**Figure 9.** Mean daily discharge: (**a**) available readings of mean daily discharge at the bridge site; (**b**) 175 years of mean daily discharge simulation using ARIMA model with p:5, d:1, q:7 and a Gaussian innovation distribution.

**Figure 10.** Simulation outcome for scour–fill depth over time (simulation starts at year 64): (**a**) one 170-year simulation; (**b**) 1064 simulations of 170 years.

The maximum scour depth value is obtained from each one of the 1064 simulations, and the best fit distribution is obtained as per 2.4 of this study. Once the probability function is found, the survival function is obtained as follows:

$$P[y\_s(te) > y \mid t] = 1 - \int\_{-\infty}^{t} \exp\left[-\left(1 + \xi \frac{t - \mu}{\sigma}\right)^{-\frac{1}{\xi}}\right] dt\tag{36}$$

Figure 11a shows the survivor function of the GEV function for different time intervals, and Figure 11b shows the PDF of the maximum scour depth.

**Figure 11.** Outcome of the simulations for different time intervals of maximum scour depth: (**a**) GEV survivor function, 100-year return period scour depth is shown as a blue line; (**b**) GEV PDF.

Results show that, as expected, scour accumulates through time. Figure 11b shows the PDFs of several time intervals of the simulation; the distribution tends to zero skewness as time increases, and it has a negative skewness and positive kurtosis at lower time intervals. A scour depth of 2.67 m was computed using the HEC 18 methodology assuming a flood with a return period of 100 years, and it is used for comparison. Such a threshold has a notable probability of exceedance, as shown in Figure 11a. Multiple events of scour with low return periods can accumulate in time, and they can be larger than the design event when the bridge is not repaired. As time tends to infinity, the maximum scour depth increases because, with the approach proposed, scour depth is not limited by harder soils that are naturally found in deeper layers.

Figure 12a shows the scour survivor function of several time intervals using the SRICOS-EFA methodology [5]. Figure 12b shows a comparison between the proposed approach that takes in to account an estimation of sediment deposit, which lessens the severity of scour accumulation in time, and the SRICOS-EFA methodology.

**Figure 12.** Simulations for different time intervals of maximum scour depth: (**a**) SRICOS EFA survivor function; (**b**) comparison between the proposed methodology and SRICOS EFA [5].

#### *5.5. Discussion*

The difference in results shown in Figure 12b are due to the differences in event simulation and in the methodology used to estimate the scour depths. Some notable differences are as follows: (**a**) SRICOS-EFA does not take in to account every event for the scour accumulation (flood has to be larger than flood to produce scour), and in a similar way a shear threshold is needed to produce scour. (**b**) It relies on average daily discharge produced by an ARIMA model, which has limitations in dealing with outliers. (**c**) It does not take in to account the accretion process (thus the sharp fall in the survivor function). (**d**) The SRICOS-EFA method uses equations obtained from flume test conducted at Texas A&M University [5], which are not equal to the HEC 18 equations for scour depth, and this leads to a natural spread on the comparison.

The proposed model has the following key aspects that differ from the SRICOS-EFA: (**a**) larger scour depths in the proposed methodology are likely due to the BMA approach on the maximum annual flow, and this approach leads to higher discharge and, thus, higher scour depths; (**b**) the addition of the accretion process leads to a heavy tail in the GEV distribution and a smoother fall in the survivor function; (**c**) in the proposed approach the hydraulic radius changes in every time step (unlike the SRICOS-EFA); (**d**)every discharge causes scour since the HEC 18 methodology does not establish a shear threshold for the scour development.

## *5.6. Scour Hazard*

Scour hazard curves provide annual exceedance probabilities of scouring events having different intensities. Since scour is time-dependent, several time interval curves are needed because such curves are developed using all the available data of the events that could trigger the scour or fill. They reflect the behavior of the phenomena in a timedependent fashion, and they can be used to evaluate both new and deteriorating bridges (Figure 13). The estimates of maximum scour are derived from simplified 1D equations. Such estimates are not exact and do not predict the future, but they provide a reasonable approximation to it.

**Figure 13.** Scour hazard curves.

Figure 13 shows different exceedance rates due to scour or scour hazard curves for different time intervals. It is noticed that the scour depth increases as the time interval increases, which means that a certain scour depth takes different values for each time interval. Scour hazard curves may be improved if the next items are included in the model. First of all, a multi-layer approach is needed to improve the estimation of the scour process, and either a nonstationary stochastic process or an autoregressive integrated moving average with additive outliers in the ARIMA-AO model is needed to calculate the discharge forecast. In addition, it is important to consider climate change, which could impact the parameters related to discharge. It is also important to verify the computed scour depth with the actual bridge scour.

## **6. Concluding Remarks**

A time-dependent scour hazard model that took into account missing data and the interaction between erosion and accretion was developed. The proposed model is straightforward, can be easily applied, and can be useful for the design and the re-design of new bridges prone to scour deterioration. In the case of existing bridge structures, the model can be used to estimate the time instant in which the structure could present an undesirable performance level caused by scour effects.

The proposed stationary stochastic model developed to study the scour and fill process shows that a single event with a high return period, such as *Tr* = 100 years, might not be adequate to represent the maximum scour that could be developed throughout the return period. Bridges subjected to continuous events of scour and fill without maintenance can acquire scour depths greater than the designed scour depth. This phenomenon makes bridges more vulnerable and partially explains why scour is the leading cause of bridge collapses in Mexico.

Scour hazard curves for different time intervals were estimated based on a reinforced concrete bridge located in one of the zones with high values of river discharge, such as Oaxaca. The scour hazard curves provide both an idea of an expected exceedance of scour associated with different time intervals and represent an important tool to estimate the reliability due to scour. Moreover, in zones with a high frequency of seismic occurrences, the scour hazard and seismic hazard can be treated from a multi-hazard point of view to estimate reliability indicators such as fragility curves, exceedance demand rates, mean annual rate of failure, or confidence factor. The above indicators lead engineers and decision-makers to re-design or repair the structural system in the case of new or existing bridges, respectively. The approach used to preprocess, characterize, and simulate events can be used for any phenomenon that has an "extreme" nature. In addition, the inclusion of missing data in the model improves the scour hazard estimation.

**Author Contributions:** Conceptualization, D.F.-V., D.T. and R.G.; methodology, D.F.-V., D.T. and R.G.; software, D.F.-V.; validation, D.F.-V., D.T. and R.G.; data curation, D.F.-V.; writing—original draft preparation, D.F.-V. and R.G.; writing—review and editing, D.F.-V., R.G. and D.T.; funding acquisition, R.G. and D.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Universidad Nacional Autónoma de México.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the manuscript.

**Acknowledgments:** The first author gratefully acknowledges economic support from Instituto de Ingeniería of Universidad Nacional Autónoma de México during his PhD studies. The second author gratefully acknowledges the support provided by the Instituto de Ingeniería of Universidad Nacional Autónoma de México.The third author appreciates the support given by both Universidad Autónoma Metropolitana and Consejo Nacional de Ciencia y Tecnología through the Ciencia Básica Project CB 2017-2018 A1-S-8700.

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

#### **References**


## *Article* **Turbulent Flow through Random Vegetation on a Rough Bed**

**Francesco Coscarella 1,\* , Nadia Penna <sup>1</sup> , Aldo Pedro Ferrante <sup>1</sup> , Paola Gualtieri <sup>2</sup> and Roberto Gaudio <sup>1</sup>**


**Abstract:** River vegetation radically modifies the flow field and turbulence characteristics. To analyze the vegetation effects on the flow, most scientific studies are based on laboratory tests or numerical simulations with vegetation stems on smooth beds. Nevertheless, in this manner, the effects of bed sediments are neglected. The aim of this paper is to experimentally investigate the effects of bed sediments in a vegetated channel and, in consideration of that, comparative experiments of velocity measures, performed with an Acoustic Doppler Velocimeter (ADV) profiler, were carried out in a laboratory flume with different uniform bed sediment sizes and the same pattern of randomly arranged emergent rigid vegetation. To better comprehend the time-averaged flow conditions, the time-averaged velocity was explored. Subsequently, the analysis was focused on the energetic characteristics of the flow field with the determination of the Turbulent Kinetic Energy (TKE) and its components, as well as of the energy spectra of the velocity components immediately downstream of a vegetation element. The results show that both the vegetation and bed roughness surface deeply affect the turbulence characteristics. Furthermore, it was revealed that the roughness influence becomes predominant as the grain size becomes larger.

**Keywords:** rigid vegetation; bed roughness; turbulent flow; Turbulent Kinetic Energy (TKE); energy spectra

## **1. Introduction**

Vegetation exerts important effects on hydraulic resistance, turbulent structures, mixing processes and sediment transport in rivers [1–5]. For this reason, a large amount of both experimental and numerical researches has been devoted to the study of the impacts of vegetation on the flow characteristics, influencing mass and momentum exchange across the river section, together with geomorphology, water quality and aquatic biodiversity (e.g., [6,7]). The flow through emergent rigid vegetation has been widely investigated, neglecting impacts introduced by natural vegetation usually observed in many fluvial ecosystems [8]. In order to understand the flow evolution in the presence of emergent and submerged vegetation, which is founded on different specific aspects of the canopy flow (mean momentum balance, turbulence budget, exchange dynamics), several works have focused on flexible vegetation (e.g., [9–11]).

Maji et al. [12] compiled a state-of-the-art study that included works on flow dynamics and interactions between flow and vegetation. Most of them aimed only at the study of the flow–vegetation interactions on smooth beds (e.g., [1,2,13–21]). Nevertheless, special interest should be devoted to works on vegetated flows with rough beds, since the interactions between fluid, vegetation and bed sediment allow for reaching better knowledge of the turbulence characteristics in real rivers, which have a crucial role in sediment transport. In fact, with respect to a smooth bed, in the case of a rough bed a pronounced velocity spike occurs near the rough surface and immediately downstream of a stem [22]. In addition, a rough bed induces a decrease in the temporal-averaged streamwise velocity and an increase in the Turbulent Kinetic Energy (TKE) [23]. In rough conditions, the streamwise velocity

**Citation:** Coscarella, F.; Penna, N.; Ferrante, A.P.; Gualtieri, P.; Gaudio, R. Turbulent Flow through Random Vegetation on a Rough Bed. *Water* **2021**, *13*, 2564. https://doi.org/ 10.3390/w13182564

Academic Editor: Mouldi Ben Meftah

Received: 30 July 2021 Accepted: 15 September 2021 Published: 17 September 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/).

also has a quasi-constant distribution in the layer where the flow is principally controlled by the vegetation [16]. Then, approaching the bed, it reduces logarithmically from the maximum constant magnitude toward zero. Moreover, the Reynolds shear stresses present very small values, practically vanishing in the region dominated by the vegetation [16,22].

Very recently, Penna et al. [24,25] studied the flow field around a rigid cylinder in three different rough bed conditions with a uniform pattern of stems that are regularly aligned. Penna et al. [24] analyzed the velocity, shear stress distributions, TKE and the energy spectra, showing that, in the region near the free surface, the flow is deeply affected by the stems. Moving toward the bed surface, the flow is influenced by both the vegetation and bed roughness effects. Penna et al. [25], using the so-called Anisotropy Invariant Maps (AIMs), investigated for the first time the turbulence anisotropy through uniform emergent rigid vegetation on rough beds. The study of the AIMs indicated that, approaching the bed surface, the combined impact of vegetation and bed roughness affects the turbulence evolution from the quasi-three-dimensional isotropy to axisymmetric anisotropy. This proved that, as the influence of the bed roughness decreases, the turbulence tends to isotropy.

Starting from the results of Penna et al. [24,25] in the case of a uniform vegetation pattern, the aim of this work is to study the impact of bed roughness and random vegetation pattern distribution on turbulence. In particular, in order to describe the flow domain, the time-averaged approaching flow velocity field, TKE, normal shear stresses and energy spectra of the velocity components were computed in an area centered on a single stem.

The paper is organized as follows: Section 2 describes the laboratory and the methodology applied for the data analysis; Section 3 illustrates and discusses the results; Section 4 reports the conclusions of the present work.

#### **2. Laboratory Experiments and Methodology**

The experimental study was conducted in a 9.6-m long, 0.485-m wide and 0.5-m deep tilting flume at the *Laboratorio "Grandi Modelli Idraulici" (GMI)*, *Università della Calabria*, Italy. In order to reduce the influence of the pump on the turbulence characteristics of the flow, a stilling tank, an uphill slipway and honeycombs (10 mm in diameter) were placed at the inlet of the channel. At the outlet, a tank equipped with a calibrated Thomson weir to measure the flow discharge *Q* and with a tailgate to regulate the water depth *h* of the flow were placed (Figure 1).

**Figure 1.** Schematic of the experimental facility (dimensions are expressed in meters).

The experiments were carried out with a flow depth *h* ≈ 0.12 m, measured 50 cm upstream of the vegetation array (i.e., in undisturbed flow condition) by a point gauge with a decimal Vernier having an accuracy of ±0.1 mm and a flow discharge *Q* equal to 19.73 L/s (measured with a Thomson weir). The approaching cross-section average flow velocity *U* = *Q*/(*Bh*) was, hence, equal to 0.30 m s−1, where *B* was the flume width. The longitudinal bottom slope of the flume, *S*, was fixed at 1.5‰ using a hydraulic jack.

The rigid vegetation was simulated with vertical, wooden and circular cylinders. The cylinder height and diameter were *hc* = 0.40 m and *d* = 0.02 m, respectively. The stems were implanted into a 1.96-m long, 0.485-m wide and 0.015-m thick Plexiglas panel fixed to the channel bottom. A total of 68 cylinders were randomly arranged in the flume (Figure 2a). All the experiments were carried out in conditions of emergent vegetation. The frontal

area per volume was *a* = *nd* = 1.4 m<sup>−</sup>1, where *n* = 71 m−<sup>2</sup> was the number of cylinders per bed area, while the solid volume fraction occupied by the canopy els per bed area, while the solid volume fraction occupied by the canopy elements was *φ* = *πad*/4 = *nπd*24 = 0.02, which is consistent with typical laboratory studies with vegetation [17,26,27]. Following Nepf [28], this vegetation distribution can be classified as dense. Focus was given to the evolution of the turbulence characteristics in a study area selected around a single stem (Figure 2b).

**Figure 2.** (**a**) Sketch of the cylinder array in the laboratory flume (dimensions are in cm); (**b**) the measurement verticals within the study area. Here, *Bs* and *Ls* are the width and length of the study area, respectively.

Three different types of bed roughness were simulated, employing very coarse sand (*d*<sup>50</sup> = 1.53 mm), fine gravel (*d*<sup>50</sup> = 6.49 mm) and coarse gravel (*d*<sup>50</sup> = 17.98 mm), respectively (Figure 3). The grain size distributions were relatively uniform, i.e., as reported by Dey and Sarkar [29], with a geometric standard deviation *σ<sup>g</sup>* = (*d*84/*d*16) 0.5 < 1.5, where *d*<sup>16</sup> and *d*<sup>84</sup> are the sediment sizes for which 16% and 84% by weight of sediment is finer, respectively. At the beginning of each run, the flume was filled in with the sediments, which were successively screeded to make the longitudinal bed slope equal to that of the flume bottom.

**Figure 3.** Sediments used in the experimental runs: (**a**) *d*<sup>50</sup> = 1.53 mm; (**b**) *d*<sup>50</sup> = 6.49 mm; (**c**) *d*<sup>50</sup> = 17.98 mm.

The experimental conditions used in this research do not reproduce a specific situation of a real river. Nevertheless, they could represent a new dataset that, for instance, may be used for the calibration of advanced numerical models. In fact, as is known, vegetation, bed roughness or man-made structures acting as an obstruction for the flow generate turbulence and affect the entire flow velocity distribution, modifying the turbulence behavior [30,31].

In Table 1, the hydraulic conditions of the experimental study are reported. Along with the aforementioned characteristics, in Table 1 the following quantities are listed: the shear velocity (*u*\*), the critical velocity for the inception of sediment motion (*Uc*), the mean water temperature (*T*) measured with the Acoustic Doppler Velocimeter (ADV) integrated thermometer (having an accuracy of ±0.1 ◦C), the water kinematic viscosity (*ν*), computed as a function of the water temperature [32], the flow Froude number *Fr* [=*U*/(*gh*) 0.5], the flow Reynolds number *Re* (=*Uh/ν*), the shear Reynolds number *Re*\* (=*u\*ε*/*ν*, where *ε* is the Nikuradse equivalent sand roughness, equal to about 2*d*50) and the Reynolds number of the vegetation stems *Red* (=*Ud*/*ν*). In accordance with Manes et al. [33] and Dey and Das [34], the shear velocity used to scale the flow statistics was determined as *u\** = (*τ\**/*ρ*) 0.5, where *τ\** is the total stress acting at the roughness tops. This can be obtained by extending linearly the distribution of the turbulent shear stress captured 50 cm upstream of the vegetation pattern (i.e., in correspondence with the undisturbed flow condition) from the region above the roughness elements to their tops. Thus, the shear velocity was evaluated at the sediment crest level as −*uw* 0.5 , where *u* and *w* are the fluctuations of the temporal velocity signal in the streamwise and vertical directions, respectively, and the symbol · indicates the time averaging operation. The critical velocity for the inception of sediment motion *Uc* was established 50 cm upstream of the vegetation array through the well-known Neill formula [35], as follows:

$$
\Delta L\_c = \sqrt{2.5 \left(\frac{h}{d\_{50}}\right)^{0.2} g \Delta d\_{50}} \tag{1}
$$

where *g* is the gravitational acceleration, Δ = (*ρ<sup>s</sup>* − *ρ*)/*ρ* is the relative submerged grain density, *ρ<sup>s</sup>* is the grain density and *ρ* is the fluid density. All the experiments were performed in clear-water condition (*U* < *Uc*), which was also verified from the direct observation of the flow.



An ADV profiler with down-looking probe, four beams (Nortek Vectrino) and an automatic movement system (the Traverse System by HR Wallingford Ltd., Oxfordshire UK) was used to capture the instantaneous velocity components (streamwise *u*, spanwise *v* and vertical *w*) with an accuracy of ±5% (assessed in previous works). The instrument sampling frequency was 100 Hz, and the duration of a single sampling was 300 s for a total number of samples of 30,000 which, as reported by [34,36,37], is adequate for determining accurate turbulence statistics. The sampling volume was a 1 mm long cylinder with a

diameter of 6 mm. The ADV receivers pointed at 50 mm below their own transmitter. Hence, the measurements were not performed near the free surface flow zone (i.e., 50 mm below the free surface). The spatial coordinates of the Traverse System had an accuracy of ±0.1 mm. Prior to the analysis of the ADV data, it was necessary to proceed with the spike detection. Firstly, the ADV raw data were prefiltered, discarding the values with correlation (COR) lower than 70% and a signal-to-noise ratio (SNR) lower than 15 dB [34]; secondly, the contaminated velocity records were cleaned using the phase-space thresholding method and each spike was replaced with a cubic polynomial through 12 points on either side of itself [38]. The de-spiking method resulted in a rejection of less than 5% of the original velocity time series.

For all the runs, 44 vertical profiles were captured, as shown in Figure 2b. The vertical spatial resolutions were 3 mm for *z* ≤ 15 mm and 5 mm above, where *z* is the vertical axis starting from the maximum crest level in the study area.

To describe the undisturbed flow 50 cm upstream of the vegetation array, the velocity vertical distribution was captured at the centerline of the laboratory flume during each run. In Figure 4a, the undisturbed profiles of the dimensionless time-averaged velocity in the streamwise direction *u*ˆ*UP* (= *uUP*/*u*∗, where *u* is the time-averaged velocity in the same direction) for the three experimental runs are reported. The vertical axis *z*ˆ was made dimensionless by dividing the elevation *z* by the local water level that, for the undisturbed distributions, was equal to the flow depth *h* reported in Table 1. As the elevation *z* increases, the streamwise velocities *u* increase; instead, near the sediment grains, in the so-called roughness sublayer, they tend to zero owing to the bed roughness (this is typical in the open-channel flow condition) [39]. In particular, for Run 3, owing to higher roughness dimension, the streamwise velocity profile tended rapidly to zero starting from elevation *z* = 0.1*h* [40]. In addition, *u*ˆ in the three runs shows different values at a fixed *z*ˆ. This happens owing to the different bed roughness conditions that lead to an increase of the shear velocity as *d*<sup>50</sup> increases and, consequently, to a decrease in *u*ˆ. The distributions of the dimensionless turbulent shear stresses *<sup>τ</sup>*ˆ*uw* (= −*uw*/*u*<sup>2</sup> <sup>∗</sup>) and of the dimensionless viscous shear stresses *τ*ˆ*<sup>ν</sup>* [= *ν*(d*u*/d*z*)/*u*<sup>2</sup> <sup>∗</sup>] along *z*ˆ are represented in Figure 4b,c, respectively. In particular, above the roughness surface, the prevalence of the Reynolds shear stresses can be noted, while the viscous shear stresses are practically negligible as *z*ˆ increases. The viscous shear stresses achieve their maximum values near the grain crests for each experimental run. Conversely, the turbulent shear stresses reach the peak above the crest level and then they reduce as the vertical distance increases.

**Figure 4.** Vertical profiles of (**a**) dimensionless time-averaged velocity, (**b**) dimensionless Reynolds shear stress and (**c**) dimensionless viscous shear stress for the undisturbed flow condition (50 cm upstream to the vegetation array) in Run 1, Run 2 and Run 3.

## **3. Results and Discussion**

## *3.1. Time-Averaged Flow*

To investigate the time-averaged flow velocity field with respect to the approaching flow velocity in the spatial flow domain (i.e., in the upstream plane section identified with letters from LE to RE in the study area), Figure 5 shows the colormaps of the dimensionless time-averaged accelerated and decelerated flow fields upstream of the investigated stem in all the runs in the plane *y*ˆ-*z*ˆ. The abscissa, represented by *y*ˆ, was made dimensionless by dividing *y* by the study area width (*Bs* = 12 cm). The time-averaged flow is accelerated if (*u*ˆ*UP* − *u*ˆ) < 0 (blue values in the colormaps) and decelerated if (*u*ˆ*UP* − *u*ˆ) > 0 (red values in the colormaps), where *u*ˆ is the dimensionless time-averaged streamwise velocity.

**Figure 5.** Contours of the dimensionless time-averaged accelerated (*u*ˆ*UP* − *u*ˆ) < 0 and decelerated (*u*ˆ*UP* − *u*ˆ) > 0 flow field in the upstream plane section (from LE to RE) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

A strong spanwise variation of *u*ˆ can be noted with respect to the approaching flow profile. Specifically, in correspondence with the cylinder (0.42*Bs* < *y* < 0.58*Bs*), in all the runs the flow is decelerated. This probably occurs owing to the incidence of both the analyzed stem and another stem immediately upstream of the former (Figure 2a). An examination of the contours also reveals a decelerated flow zone on the left side of the cylinder (0< *y* < 0.2*Bs*) in all the experiments. This happens owing to the presence of another stem immediately upstream of the study area (Figure 2a). Instead, the flow fields are accelerated both on the right side of the analyzed cylinder and, to a lesser extent, on the left one. Moving toward the bed, it is evident that the flow velocity is also influenced by the bed roughness. As *d*<sup>50</sup> decreases, this zone becomes strongly accelerated. Conversely, as *d*<sup>50</sup> increases, in the near-bed layer, the flow field results to be influenced by the bed roughness and much more by the vegetation, with a minor acceleration intensity.

A similar behavior can be appreciated in the downstream flow domain in Figure 6 (i.e., in the downstream plane section identified with letters from LV to RV in the study area). In particular, immediately downstream of the cylinder (0.42*Bs* < *y* < 0.58*Bs*), the flow is decelerated in all the runs. Conversely, Run 2 shows an accelerated flow both on the right and on the left sides of the cylinder, while Run 1 and Run 3 show an accelerated flow mostly on the right side of the cylinder. A sensible difference is clear in the values of *u*ˆ*UP* − *u*ˆ immediately upstream (vertical A; *y* = 0.5*Bs*) and downstream (vertical W; *y* = 0.5*Bs*) of the cylinder. Specifically, the time-averaged streamwise velocity *u* is reduced by about 30% and 50% with respect to the time-averaged streamwise velocity of the undisturbed profile, *uUP*, at the verticals A and W, respectively, in all the runs. It is possible to assure that a major deceleration is obtained beyond a cylinder. In fact, as observed in the downstream section at *y* = 0.5*Bs* (vertical W) with respect to the upstream section at *y* = 0.5 *Bs* (vertical A), the velocity reduction is practically due to the vicinity of the studied stem (6 cm upstream). On the contrary, a minor deceleration is visible at vertical A, although it is between two cylinders (an upstream stem at 9 cm and a downstream stem at 6 cm).

**Figure 6.** Contours of the dimensionless time-averaged accelerated (*u*ˆ*UP* − *u*ˆ) < 0 and decelerated (*u*ˆ*UP* − *u*ˆ) > 0 flow field in the downstream plane section (from LV to RV) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

In order to investigate the longitudinal evolution of the flow field, the contours of the dimensionless time-averaged accelerated and decelerated flow field are analyzed. For the sake of simplicity, they are reported in Figure 7 only for the extreme right vertical plane identified with letters from RE to RV in the plane *x*ˆ-*z*ˆ. The abscissa *x*ˆ was made dimensionless by dividing *x* by the study area width (*Ls* = 12 cm). In Figure 7, it is evident that the presence of vegetation has a very visible effect on the flow field: the flow is accelerated in all the runs with higher streamwise velocities than in the undisturbed flow profiles along the whole water depth. At each measurement location, the vegetation causes the velocity profile to maintain a constant value (for *z* > *hl*). Moving toward the bed (*z* < *hl*), the influence of vegetation decreases, and the flow field becomes more accelerated, owing to the presence of the bed.

**Figure 7.** Contours of the dimensionless time-averaged accelerated (*u*ˆ*UP* − *u*ˆ) < 0 and decelerated (*u*ˆ*UP* − *u*ˆ) > 0 flow field in the longitudinal plane section (from RE to RV) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

## *3.2. TKE and Normal Stresses*

The computation of the spatial distributions of the TKE is very significant in the assessment of the energetic process in open-channel flows. Since the time-averaged spatial fluctuations influence mechanical dispersion [41] and, in turn, this latter may be influenced by the vegetation stems and the bed roughness, the topic is of considerable interest. In fact, the presence of vegetation adds a further turbulence production in the wakes of the plant elements [42]. The TKE is defined as half the sum of the variances of the velocity components:

$$\text{TKE} = \frac{1}{2} \left[ \overline{\left(u'\right)^2} + \overline{\left(v'\right)^2} + \overline{\left(w'\right)^2} \right]. \tag{2}$$

where *v* is the velocity fluctuation of the spanwise velocity *v*.

The colormaps of the dimensionless TKE (i.e., TKE divided by *u*<sup>2</sup> <sup>∗</sup>) on the transversal upstream plane section (identified with the letters from LE to RE in the spanwise direction) are shown in Figure 8 for each experimental run. The highest values of the TKE are located in front of the cylinder and at the vertical LE (*y* = 0).

**Figure 8.** Contours of the dimensionless TKE in the upstream plane section (from LE to RE) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

To analyze this behavior, in Figure 9 we consider, for all the runs, the dimensionless normal stresses in the streamwise, spanwise and vertical directions, *uu*/*u*<sup>2</sup> <sup>∗</sup>, *<sup>v</sup>v*/*u*<sup>2</sup> <sup>∗</sup> and *ww*/*u*<sup>2</sup> <sup>∗</sup>, respectively, which are essentially the three addends in Equation (2). It is possible to observe that the high value on the left part of TKE contours (*y* = 0 in Figure 8a) is practically ascribable to the streamwise and spanwise effects (Figure 9) owing to the presence upstream of the studied stem, both of an empty zone without vegetation, which influences *uu* , and of a wake vortex, that affects *vv* . In contrast, the high magnitude of the TKE immediately behind the vegetation element (0.42*Bs* < *y* < 0.58*Bs*) is mostly due to the spanwise fluctuations, as a consequence of the circumvention of the obstacle.

**Figure 9.** Contours of the dimensionless normal stresses (**a**) *uu*/*u*<sup>2</sup> <sup>∗</sup>, (**b**) *<sup>v</sup>v*/*u*<sup>2</sup> <sup>∗</sup> and (**c**) *<sup>w</sup>w*/*u*<sup>2</sup> <sup>∗</sup> in the upstream plane section (from LE to RE) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

Conversely, downstream of the studied vegetation element, the dimensionless TKE shows higher magnitudes immediately beyond the stem, whereas no high TKE value is detectable at *y =* 0 (Figure 10).

**Figure 10.** Contours of the dimensionless TKE in the downstream plane section (from LV to RV) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

Analogously to the upstream flow field, from the normal shear stress distributions shown in Figure 11, it is possible to evaluate the contributions of the TKE patterns. Specifically, the higher kinetic energy values are influenced by both the streamwise, *uu* , and spanwise, *vv* , normal stresses, which increase owing to the presence of the von Kármán wake vortex. This behavior, although with different magnitudes, is displayed in all the runs and, consequently, is clearly a vegetation effect. Furthermore, from a comparison between the TKE and the normal stress lateral distributions (Figures 8 and 10, and Figures 9 and 11, respectively), it is evident a predominant effect of the vegetation element in the downstream plane section and a contribution of the vertical normal stresses in the vertical W (downstream of the cylinder) greater than in the vertical A (upstream of the cylinder).

**Figure 11.** Contours of the dimensionless normal stresses (**a**) *uu*/*u*<sup>2</sup> <sup>∗</sup>, (**b**) *<sup>v</sup>v*/*u*<sup>2</sup> <sup>∗</sup> and (**c**) *<sup>w</sup>w*/*u*<sup>2</sup> ∗ in the downstream plane section (from LV to RV) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

The distributions of the dimensionless TKE in the longitudinal extreme right vertical plane (identified with letters from RE to RV) are illustrated in Figure 12. High TKE magnitudes are observed in the near-bed flow zone, where the bed roughness surface causes higher fluctuations of the velocity components. However, the TKE reduces progressively

moving upwards from *z* > 0, owing to the inhibition of *u* , *v* , and *w* . The TKE contours are slightly spatially nonuniform for all the experiments. This is probably due to the random vegetation array and the spatial irregular bed sediment.

**Figure 12.** Contours of the dimensionless TKE in the longitudinal plane section (from RE to RV) for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3. The black dashed lines show the cylinder position.

## *3.3. Energy Spectra*

In order to further analyze the turbulent characteristics, the measured velocity data were explored through the energy spectra of the velocity fluctuations. In particular, the energy spectra are shown in Figure 13 for the vertical W, i.e., immediately downstream of the studied cylinder for all the runs and three different elevations (*z* = 0, *z* = 0.2 *h* and *z* = 0.4 *h*), as a function of the Strouhal number of the cylinder (*St* = *fd*/*u*, where *f* is the frequency with a resolution equal to *Fs*/*N*, and *N* is the number of samples equal to 30,000 for an acquisition time of 300 s). The energy spectra were determined by employing the discrete fast Fourier transform of the autocorrelation function.

**Figure 13.** Energy spectra for (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3 of streamwise (blue lines), spanwise (red lines), and vertical (green lines) velocity fluctuations at three different levels *z*ˆ of the vertical W.

From a comparison of the spectra, the spanwise velocity revealed the presence of large-scale coherent structures, evident as a peak located in the energy-containing range of

the energy spectra (Figure 13). Specifically, these peaks are present in all the runs, with the exception of two spectra of Run 3, where it is evident that the peak recedes as *z* decreases. This may be due to a roughness effect, which in Run 3 is explicated with a medium sediment diameter equal to 17.98 mm, i.e., comparable with the stem diameter, equal to 20 mm. Moreover, the peaks were observed at the same Strouhal number in the *u*'*v*' crossspectra (not shown here for the sake of brevity), demonstrating that the identified coherent structures were responsible across the vegetation for the lateral momentum transport [5].

#### **4. Conclusions**

In the present work, an experimental study was carried out to investigate the impact of different uniform bed roughness on the flow characteristics through randomly arranged emergent rigid vegetation. In particular, focus was given to the evolution of the turbulence characteristics in a study area selected around a single stem. The principal results are summarized below.

In the sections, respectively, upstream and downstream of the generic stem, the dimensionless time-averaged accelerated and decelerated flow field is deeply affected by the vegetation and, to a lesser extent, by the bed roughness. In particular, in the flow field downstream of the studied cylinder it is possible to point out these zones: a decelerated one, in correspondence to the stem (0.42*Bs* < *y* < 0.58*Bs*), another decelerated zone, on the left side of the stem (0 < *y* < 0.2*Bs*, only for Runs 1 and 3) and an accelerated one elsewhere.

From a comparison perspective, it is possible to evaluate a velocity reduction of about 50% in the downstream plane immediately behind the analyzed stem. Conversely, downstream of the cylinder, on the left side (0< *y* < 0.2*Bs*), the decelerated zones fade about 30%. This is clearly attributable to the random distribution of the vegetation pattern. Instead, in the longitudinal plane the flow is decelerated in all the runs with higher streamwise velocities than in the undisturbed flow profiles along the whole water depth.

The analysis of the TKE distribution clearly shows the effects of the vegetation, with high magnitudes immediately upstream and downstream of the investigated stem. This indicates that the velocity oscillations get excited by the cylinders, producing an increased turbulence intensity in the proximity of the stem. Moving toward the free zone (i.e., without vegetation), this influence vanishes, causing a decrease in the TKE value. In addition, the behavior observed in TKE colormaps was examined by comparing the streamwise, spanwise and vertical normal stress contours. The results revealed a strong influence of the *u*' and *v*' fluctuations on the energy distribution and highlighted the influence of the so-called von Kármán vortices. In addition, from the analysis of the TKE distributions, an effect of the cylinder greater in the downstream plane section than in the upstream one is manifested, with an increased mean value of 25% at the abscissa *y* = 0.5 *Bs*. The longitudinal TKE distribution revealed high values in the near-bed flow zone. This suggests that the velocity oscillations get excited by the rough bed, producing an increase of the turbulence level in the vicinity of the sediments. Mowing toward the free surface this effect disappears, inducing a decrease in the TKE.

The evaluation of the energy spectra of the velocity fluctuations showed a clear influence of both vegetation and roughness. In particular, the spanwise velocity component revealed energetic peaks that indicate the presence of large-scale coherent structures due to the von Kármán wake vortex along the flow depth in the case of lower roughness (Run 1 and 2) and only in the upper part of the water depth in the case of higher roughness (Run 3). In fact, it is interesting to point out that, when the median sediment diameter is comparable with the stem diameter (i.e., Run 3) going towards the bed, the energetic peak is lowered. This involves a consideration of the nature of coherent structures of turbulence, which are significantly influenced by the characteristic scales of the flow conditions, such as vegetation diameter, water depth and roughness size.

Further study is necessary to describe more deeply the turbulence structures in the presence of both vegetation and sediments. In fact, the results referring to a single generic

cylinder are representative only for such a case and cannot be generalized, owing to the lack of further data relating to different vegetation elements. For this purpose, in future works we intend to perform other laboratory experiments, intensifying the measurements and adopting advanced techniques, such as the Particle Image Velocimetry (PIV), in order to better characterize the two-dimensional (2D) turbulence structures.

**Author Contributions:** Conceptualization, F.C., N.P., A.P.F., R.G.; methodology, F.C., N.P., P.G., R.G.; formal analysis, F.C., N.P.; data curation, F.C.; writing—original draft preparation, F.C.; writing review and editing, F.C., N.P., A.P.F., P.G., R.G.; supervision, R.G.; funding acquisition, R.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the "*SILA*-*Sistema Integrato di Laboratori per l'Ambiente*-An Integrated System of Laboratories for the Environment" PONa3\_00341.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. This data can be found here: http://www.ingegneriacivile.unical.it/lgmi/turbulent-flow-through-randomvegetation-on-a-rough-bed/, accessed on 14 September 2021.

**Acknowledgments:** The authors would like to thank Erik Calabrese Violetta for his valuable work during the performance of the experimental runs.

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

## **References**

