**Man-Induced Discrete Freshwater Discharge and Changes in Flow Structure and Bottom Turbulence in Altered Yeongsan Estuary, Korea**

**KiRyong Kang <sup>1</sup> and Guan-hong Lee 2,\***


Received: 3 June 2020; Accepted: 1 July 2020; Published: 5 July 2020

**Abstract:** Flow measurements were performed in the altered Yeongsan estuary, Korea, in August 2011, to investigate changes in flow structure in the water column and turbulence characteristics very close to the bed. Comparison between the bottom turbulent kinetic energy (TKE) and suspended sediment concentration (SSC) was conducted to examine how discrete freshwater discharge affects the bottom sediment concentration. The discrete freshwater discharge due to the gate opening of the Yeongsan estuarine dam induced a strong two-layer circulation: an offshore-flowing surface layer and a landward-flowing bottom layer. The fine flow structure from the bed to 0.35 m above the bottom (mab hereafter) exhibited an upside-down-bell-shaped profile for which current speed was nearly uniform above 0.1 mab, with the magnitude of the horizontal and vertical flow speeds reaching 0.1 and 0.01 m/s, respectively. The bottom turbulence responded to the freshwater discharge at the surface layer and the maximum magnitude of the Reynolds stress reached up to 2 <sup>×</sup> <sup>10</sup>−<sup>4</sup> m2/s2 during the discharged period, which coincided with increased SSC in the bottom boundary layer. These results indicate that the surface freshwater discharge due to opening of the estuarine dam gate increases the SSC by the discharge-induced intensification of the turbulent flow in the bottom boundary layer.

**Keywords:** Yeongsan estuary; freshwater discharge; two-layer circulation; Reynolds stress; bottom turbulence; suspended sediment concentration

#### **1. Introduction**

In estuaries, water currents are mainly driven by tide, wind, and freshwater discharge [1]. Other factors affecting the flow structure include density stratification caused by temperature and salinity, seabed roughness due to the bottom surficial texture and bedforms, and also artificial alterations such as dam construction. Water currents play an important role in the transportation and distribution of suspended sediment which can change the topography and morphological shape of an estuary. Especially when the bottom surface is composed of silt or mud, resuspension of sediment can be very dependent on the turbulent flow activity above the bed. Turbulent flow near the bed is a well-known factor which generates sediment-related processes such as erosion, dispersion, and transportation in the bottom layer [2–4].

Artificial alteration of an estuarine environment, such as dam construction to block the saltwater intrusion and the regulation of the freshwater discharge, can modify the physical environment such as the tidal range and circulation structure [5,6]. One factor which can induce a sudden response from an estuary is forced freshwater discharge because rapid currents can be formed in the surface with a velocity difference between the surface and bottom layers. When any such artificial change occurs such as discontinuous freshwater input, an estuary naturally will adjust to these changes. An example of an estuary which has experienced such artificial change is the Yeongsan River estuary (Figure 1), located in the southwestern portion of the Korean peninsula. Yeongsan River is one of the four major rivers in Korea. It has a drainage area of 3455 km<sup>2</sup> and a length of 129.5 km. The Yeongsan estuary has a series of offshore islands and extensive salt marsh complexes, with an average water depth of 20 m. The main channel was blocked by a dam in February 1981, and now it behaves as a semi-enclosed bay [7]. The tide is macrotidal with a semi-diurnal tidal range of 4.5 m, and 70% of the annual rainfall of 840 mm occurs from June to August.

**Figure 1.** Study area and location of current meters in Yeongsan estuary, Korea. Black dots mark the location of two current meters placed in August 2011.

According to several studies, the physical environmental conditions have changed after the dam construction. Kang [6] showed that the tidal range increased to 60 and 43 cm for extreme high and low tide, respectively, while the tidal velocity decreased with an ebb tidal dominance causing a change in the sediment transport mechanisms. During the summer season when the gate opens frequently, the freshwater discharge has become an important factor to change the current, temperature, and salinity distribution because of sudden and forced discharge through the surface layer [8–12]. In terms of flow system change, for example, Cho et al. [11] suggested that there exist four layers under low discharge conditions during the summer season, showing seaward flow in the surface and middle layer, and landward flow in the bottom layer and between the surface and middle layers. With freshwater discharge due to the gate opening, a two-layer circulation is formed with strong stratification between the offshore-flowing surface layer and the landward-flowing bottom layer. So, the freshwater discharge can affect directly the estuarine circulation system. In addition, through 210Pb and 7Be radioisotope geochronology in the Yeongsan estuary, Williams et al. [13] showed that high sedimentation rates up to 9 cm/year occur in the estuary, and the sediment deposition primarily occurs during episodic events corresponding to high discharge. However, the flow structure and the associated sediment suspension near the bed during freshwater discharge are still not well understood.

In the Yeongsan estuary, a sudden release of freshwater in the surface layer due to opening of the estuarine dam gate can cause rapid intensification of the seaward flow. It is expected, then, that the flow of the lower layer should respond to the sudden and strong surface flow. In such a case, the bottom flow response could be coupled to change in the velocity shear, and this in turn could be a cause for sediment resuspension from the bed. To our knowledge, this is the first field observation to report the characteristics of flow, turbulence, and suspended sediment in the bottom boundary layer in response to discrete freshwater discharge. Somewhat related to this topic includes those of entrainment and mixing in a turbulent jet [14,15], jet scour [16], and sediment discharge by jet-induced flow [17–19]. Therefore, the main objectives in this study are to investigate how the freshwater discharge changes the structure of the bottom turbulent flow properties like the Reynolds stress and the turbulent kinetic energy (TKE), and to elucidate how the bottom turbulent flow interacts with the suspended sediment near the seabed by using observed water current and suspended sediment concentration data to capture the turbulent structure and suspended sediment concentration in the bottom boundary layer during a freshwater discharge. The observation and data processing scheme are described in Section 2, and the analysis results about the evolution of the mean and turbulent flow and the bottom suspended sediment concentrations are shown in Section 3. The role of man-induced freshwater discharge on the flow structure and suspended sediment concentration fluctuations above the seabed is discussed in Section 4, and finally, a very brief conclusion and meaning of this study are shown in Section 5.

#### **2. Observation and Data Processing**

The experiment campaign was designed specifically to observe the flows in both the bottom and upper layers during the period of discrete freshwater discharge in the Yeongsan estuary, Korea, during August 2011 (Figure 1). Two current meters, an ADCP (RDI 1200 kHz, Teledyne, Poway, CA) and an AquaDoppHR (Model: AQP 9116, Nortek, Boston, MA), were moored near the Yeongsan estuarine dam in the inner estuary (Figure 1 and Table 1). The ADCP was mounted on the bed with an up-looking orientation to measure the flow profile in the water column, and the AquaDoppHR was bottom-mounted at about 1 m above bottom (mab hereafter) with down-looking orientation for near-bed turbulence measurements. The current profiles obtained from the ADCP with high percent good values of 85 and above were averaged by burst with a 30-min interval for the flow structure in the water column (Table 1).


**Table 1.** Measurement scheme for the current meters in Yeongsan estuary.

The ADCP current data were not rotated into the along- and cross-channel directions because the main direction of the freshwater discharge was east–west. The observed current data were separated into tidal and residual flows based on the tidal harmonic analysis developed by Foreman [20] and implemented into MATLAB as T\_TIDE [21]. From the residual flow, it was possible to check how the upper layer flow responded to the sudden freshwater discharge from the gate opening of the dam.

The AquaDoppHR was programmed for burst sampling to observe the fine structure of the mean and turbulent flow characteristics very near the bed. The burst sampling lasted for 8.5 min in 30-min burst intervals the same as the ADCP, however, the AquaDoppHR sampling frequency was higher at 4 Hz. The AquaDoppHR profiled from the bed to 0.35 mab with 5 cm bin size spatial resolution. The bed was detected from the trend of the acoustic signal strength along each beam, in which the signal strength generally shows maximum value at the bed. In this case, since three beams were looking at the bottom, the maximum value was located in different bins of each beam such as the 7th, 8th, and 9th bins. A possible reason for different bin numbers showing the maximum signal strength could be the bed status. Since the bed surface that the AquaDoppHR was looking down upon was not perfectly flat, and because it is hard to know the true bed status at the midpoint between the three beams, we selected the middle value of bins as the bed level.

The near-bed velocity data obtained from the AquaDoppHR were first despiked using an averaging method after visual inspection of all bursts. Then, the mean, variance, and covariance values of the different flow components of all bins for each burst were calculated to estimate the turbulent flow characteristics such as the Reynolds stress components and the turbulent kinetic energy (TKE). The Reynolds stress and the TKE were calculated by using the following two equations (Equations (1) and (2)):

$$R\_{13} = -\preccurlyeq \iota' \iota v' \succcurlyeq R\_{23} = -\preccurlyeq v' \circ \text{w'} \tag{1}$$

$$\text{TEE} = (\text{ +\text{ +\text{)/2 \tag{2}$$

where *u*- , *v*- , and *w* are the turbulent components of the east–west (*u*), north–south (*v*), and up–down (*w*) components of a velocity vector, u = (*u*, *v*, *w*), for each burst sampling period, respectively. U and V are the mean flow components of each burst for horizontal flow, for example, u = U + *u* and v = V + *v*- , and the bracket (< >) means time averaging for a burst period. Finally, the suspended sediment concentration (SSC) was calculated by conversion of the acoustic backscatter signal of the four beams of the moored ADCP. A more detailed procedure for this calculation can be found in Park and Lee [22].

#### **3. Results**

#### *3.1. Freshwater Discharge and Flow Structure Change*

One of the key factors to affect the estuarine flow in an altered estuary is the freshwater discharge since it makes both a forced seaward current in the surface layer and a density difference between the upper and lower layers. Moreover, the strength of the surface current as well as the density gradient depends on the amount of freshwater discharge. In this section, the response of the surface water after the freshwater discharge is depicted using the ADCP data. A time series of discharge is compared with the current structure and then with the change in residual flow (component with tides removed, using 17 tidal constituents) in the water column.

Figure 2 shows the discrete freshwater discharge and mean current velocity for U and V, during 5–28 August 2011. According to typical gate operation procedure, the gates were opened only during the low tide to prevent saltwater intrusion. On 11 August, over 6 <sup>×</sup> 107 tons of freshwater were discharged. After the discharge, the westward surface current rapidly increased to greater than 0.5 m/s and affected at least 5 m below the surface. At the same time, an opposite flow to the east formed in the lower layer. The sudden release of freshwater also increased the north–south component of the flow in a short period of time. When there was no freshwater discharge, the current speed was less than 0.25 m/s and the flow structure simply repeated the flood and ebb states such that eastward flow occurred during flood, westward flow occurred during ebb, and almost zero flow occurred during both high and low slack tides.

The tide-removed residual flow showed in detail how the mean flow structure could be changed by the amount of freshwater discharge (Figure 3). During freshwater discharge, the residual flow showed a two-layer system as expected wherein the upper layer moves seaward and the lower layer flows landward. The flow was seaward from the surface to 0.6 z/H (i.e., z = 9.6 m for H = 16 m) on 10 August 2011, and the depth of the seaward flow decreased as the freshwater discharge decreased. When the amount of discharged water increased, the lower layer also showed a stronger response of landward flow just after the discharge. Of course, the peak speed of the residual u-component appeared in the surface. During no or very weak discharge periods (20–22 August in Figure 3), the residual flow was also weak (almost less than 0.01 m/s) and the vertical structure was hard to specify because of the multilayered structure. This stagnant vertical structure changed with a rapid increase in the speed of the residual current when freshwater discharges occurred.

**Figure 2.** Time series of the freshwater discharge and water flow. (**a**) Freshwater discharge and water depth, (**b**) u-component of current velocity, (**c**) v-component of current velocity. In (**b**,**c**), negative values indicate west- and south-directed velocities, respectively.

The averaged residual flow structure for the whole observation period was a typical two-layer flow system that should be induced by vertical gravitational circulation in an estuary [1,23]. The mechanism of this structure is known very well, and there are two forcing factors: the barotropic and baroclinic forcing. The sea level difference between the upper and lower estuary acts as the barotropic forcing, making seaward flow in the upper layer, while the density difference forms the baroclinic forcing, making landward flow in the lower layer. The v-component showed a single-layer structure with very weak flow to the north, and the speed was less than 2 cm/s. This could have been due to the position of the dam gate which is located at the southern part of the dam (Figure 1). The two-layer structure during the entire period was an important aspect of the overall vertical structure of the flow since there was a stagnant, multilayer flow structure during the weak or no freshwater discharge periods (Figure 3d,e).

The strength of the upper- and lower-layer flow depended on the amount of the freshwater discharge. Figure 3 shows the vertical structure change of the residual flow with variation in the amount of discharge, for example, relatively high, low, and no freshwater cases. During the high freshwater discharge period on 8–11 August 2011, the maximum speed of the seaward residual flow was about 0.1 m/s in the surface layer, and a landward flow existed below 0.7 z/H with a quarter of the speed of the surface layer. As the amount of freshwater discharge decreased, the thickness of landward flow in the subsurface layer was increased from 0.7 z/H to 0.8 z/H, and the speed of the surface layer rapidly decreased. This represents the seawater response inside the estuary after blocking the freshwater release. When the barotropic forcing completely disappears due to no freshwater discharge, the water column has a multilayer flow structure, for example, a three- or four-layer structure, indicating that there were many local flows without a dominant flow pattern (Figure 3e). This indicated that during the no freshwater discharge period in this estuary, the residual flow pattern was indistinct and exhibited weak and vertically variable horizontal currents in a stagnant environment. Therefore, the freshwater discharge plays the major role in intensifying the two-layer system, including the residual circulation, in this estuary.

**Figure 3.** Change of vertical structure of the residual (tide-removed) flow: (**a**) east–west component, (**b**) north–south component, and mean flow structure of each period such as (**c**) high (i.e., Q <sup>&</sup>gt; <sup>4</sup> <sup>×</sup> <sup>10</sup><sup>7</sup> ton on 8–11 August), (**d**) low (i.e., Q < 1 <sup>×</sup> 10<sup>7</sup> ton on 12–15 August), and (**e**) no freshwater discharge (20–22 August). Two profiles on the right side of (**a**,**b**) are the residual over the entire period. The negative values on (**a**,**b**) indicate west- and south-directed residual velocities, respectively.

#### *3.2. Bottom Boundary Layer Flow and Reynolds Stress*

In the previous section, we examined the impact of freshwater discharge on the flow of the water column above the bottom boundary layer (BBL). In this section, the response of the bottom boundary layer flow to the rapid seaward flow in the upper layer is examined. Data from the AquaDoppHR, moored in a downward-looking orientation on a tripod during 24–28 August 2011, were used to show the detailed structure of the horizontal and vertical components of the BBL flow, the Reynolds stress components, and the turbulent kinetic energy.

Figure 4 shows the time variation of sea level, the amount of freshwater discharge, and the horizontal and vertical flows at the bottom. The magnitude of flow was about 0.1 m/s in the horizontal component and about 0.01 m/s in the vertical component, and the ebb and flood patterns were not clear since there were several irregular changes in the flow direction even during the flood or ebb tides. The vertical structure of flow indicated that the current magnitude of each profile was similar at 0.15–0.35 mab (Figure 4). However, it began to decrease from 0.15 mab and became almost zero near the bed due to the bottom friction. This phenomenon became clearer as the speed increased. The bottom seaward flow during the low tide was evidence of the fact that the freshwater release can directly affect the bottom flow strength and reverse the flow direction. When freshwater was discharged in the surface layer, the bottom seaward flow was intensified for a short time and then a landward flow appeared (Figure 4b). During 24–28 August 2011, the dam gates opened four times during the low tide: two openings on 24–25 August discharged about 6 <sup>×</sup> 106 tons of water and the other two openings released over 12 <sup>×</sup> 106 tons of water (Figure 4a). As the amount of freshwater discharge increased (i.e., on 26 August and 27 August), the change in the flow direction and strength was sharper than the other two cases. One possible reason that the reversed (landward) flow occurred in such a short time could be a result of water mass conservation. When freshwater moved rapidly seaward due to the dam gate opening, it could have produced a return flow in the bottom layer to make the water mass balance inside the estuarine area. The response of the north–south and up–down flow showed a spike-like change just after the freshwater discharge (Figure 4c,d). It is noted, on the other hand, that the flow speed was very weak or almost zero during the no freshwater discharge periods on 25 August and 27 August. It is also interesting to mention the fluctuation of the vertical flow since it displayed directional change (positive to negative) as the freshwater rapidly moved seaward in the upper layer even though the speed was one order of magnitude less than that of the horizontal flow (Figure 4d). The speed of the vertical flow was about 1 cm/s. This fluctuation of the vertical flow occurred during the low tide and was intensified when the freshwater was discharged. Therefore, the freshwater release by opening of the dam gate caused a rapid and strong seaward surface flow, and its impact could reach the bottom and influence the bottom flow structure such as with sudden direction changes in both the horizontal and vertical flow.

At the bottom boundary layer, the increase of near-bed flow can affect the velocity shear and thus intensify the turbulent flow activity. Figure 5 shows the variation of the Reynolds stresses and turbulent kinetic energy in the bottom boundary layer. The vertical structure is similar to the shape of an upside-down bell. The stress terms, −<*u*- *w*- >, showed a symmetric structure except for a few cases where the values became strongly negative (Figure 5a), and −<*v'w'*> displayed mostly positive values (Figure 5b). This was related to the fluctuation pattern of the flow. During the freshwater discharge, as noted above, the seaward flow was intensified for a short period of time and then the flow direction was reversed, making −<*v*- *w*- > symmetric. However, a mostly southward v-component flow was responsible for the positive values of −<*v*- *w*- >. Thus, the flow change due to the surface freshwater discharge was linked to the bottom stress intensification and turbulent flow change.

Figure 5c displays the time variation of turbulent kinetic energy. The vertical distribution of TKE is similar to the upside-down half-bell shape and the magnitude rapidly increased from the bed to 0.1–0.2 mab and then remained relatively constant above 0.2 mab. With respect to time, the TKE also increased during the freshwater discharge during low tide, and was proportional to the amount of freshwater discharge. The maximum value during this study period appeared during the late low tide on 27 August.

#### *3.3. Suspended Sediment Concentration and TKE*

The major sources of suspended sediment in estuarine environments in general are from the upstream and the sea, as well as bed erosion by flow–bed interaction. In this section, the temporal change in SSC profiles after the freshwater discharges is investigated and the relationship between the turbulent flow activity and SSC in the bottom boundary layer is examined.

Figure 6 illustrates the time variation of SSC profiles during 24–28 August 2011. There were clearly high SSCs found in the surface layer. During the second low tide of each day that freshwater discharge occurred, higher SSCs were observed. It appeared to be due to direct input of suspended sediment discharged from the upstream. As no peaks of SSC were observed when there was no freshwater

discharge, for example, during the first low tide of each day, it is evident that the freshwater discharge due to the opening of the dam gate was the main source of the suspended sediment in the surface layer. Another interesting aspect is that high SSC occurred after the freshwater discharge during the flood tide. It is likely that offshore-advected high SSCs during the low tide returned back into the estuary during the flood tide, resulting in high SSCs.

**Figure 4.** Time variation of the flow structure of the bottom boundary layer with freshwater discharge: (**a**) sea level and freshwater discharge, (**b**) east–west component, (**c**) north–south component, and (**d**) vertical component. The solid line in each panel indicates the sea surface level variation pattern and the scale was adjusted.

Below the surface layer, for example, around 10–12 mab, the SSC was mostly less than 0.01 g/L except during the second flood tide of each day. However, the SSC in the bottom layer suddenly increased and appeared to propagate into the higher layers whenever the surface freshwater was discharged. For example, SSC greater than 0.02 g/L (light green in Figure 6) reached above 6 mab during 25–27 August. In the previous section, we saw the rapid increase of bottom flow speed and stress-related terms (see Figures 4 and 5), and this kind of change in the physical factors could induce erosion of bottom sediments if the bed of this area consisted of very fine sediment such as silt, clay, or mud. According to Kim et al. [24] and Williams et al. [13], the sediments of this estuary consist mainly of silt–clay mixtures in which silt is distributed in the shallow areas and clay exists in the relatively deep areas of the central estuary. Certainly, the bed shear stress observed during the freshwater discharge (see Figure 5) can resuspend the fines high into the water column (Figure 6). Furthermore, Bang et al. [25] simulated sediment transport by using a numerical model, and showed

that the discharge-induced estuarine circulation could cause overall silt-size sediment deposition and sustained suspension of clay sediment.

The time series of the TKE and SSC in the bottom boundary layer is illustrated in Figure 7. It is clear that the SSC rapidly increased and then gradually decreased until the next event occurred. The TKE exhibited a peak during each freshwater discharge and was proportional to the amount of freshwater discharge. Likewise, the SSC matched well with the amount of freshwater discharge. This kind of synchronization was repeated during every freshwater discharge. Thus, it could be generalized from this relationship that when the turbulent flow was intensified by the strong seaward surface flow, a relatively large amount of sediment could be resuspended from the bed.

**Figure 5.** Reynolds stress and turbulent kinetic energy profiles near bed: (**a**) east–west component, (**b**) north–south component, and (**c**) turbulent kinetic energy (TKE). Black and white solid lines in the left panel indicate the sea level variation.

**Figure 6.** Variation of suspended sediment concentration during 24–28 August 2011.

**Figure 7.** Comparison between the turbulent kinetic energy (TKE) and suspended sediment concentration (SSC) changes in the bottom boundary layer.

#### **4. Discussion**

In estuaries altered by an estuarine dam, the opening of the dam gates results in a sudden release of a significant amount of freshwater to the river mouth area, resulting in changes to the physical and environmental conditions such as the formation of a strong surface current and sediment suspension. Measurement of the water currents during August 2011 in the Yeongsan estuary was carried out to

investigate the changes in the vertical flow structure, the bottom turbulent flow, and the relationship between the bottom turbulent flow activity and the variation of suspended sediment concentration during the freshwater discharge.

The freshwater discharge due to dam gate opening could significantly affect the vertical structure of flow, inducing velocity shear and turbulent flow activity in the bottom boundary layer. During the period of no or weak discharge, the flow was mostly tidal motion and the vertical structure showed multiple layers. Cho et al. [11] showed that the spreading of warm freshwater discharged over the pre-existing surface water and the intrusion of warm saline water along the bottom from the open sea produced a multilayer structure, and the multilayer structure remained throughout the summer due to strong stratification and a weak tidal current in the Yeongsan estuary. During the freshwater discharge, however, a two-layer system developed with seaward flow in the surface and landward flow in the lower layer [7,26]. The thickness of the seaward flow was affected by the amount of freshwater, and the thickness of landward flow in the subsurface layer increased as the freshwater amount decreased. The magnitude of the residual flow was about 10 cm/s in the surface during the discharge periods, and became stagnant with weak, vertically variable horizontal flows during no discharge periods.

The bottom flow near the bed did not follow the general flood and ebb cycle in tidal motion. The seaward and landward flows were, however, formed after freshwater discharge, for example, on 25–27 August. This appears as strong evidence that the freshwater release due to opening of the dam gate affected the bottom flow strength. One more interesting aspect is that the reversed (landward) flow followed after the initial seaward flow. A two-layer circulation system was suggested wherein, following the water mass conservation principle, a landward flow could happen in order to recover the water overtransported to the sea by the suddenly intensified seaward flow which could cause additional seaward movement. Considering the concomitant vertical flow pattern of up and down, the upward flow above the bed was first shown when freshwater was released and then downward flow happened with a maximum magnitude of 0.01 m/s.

An impact of surface freshwater discharge to the bottom boundary layer was an intensification of the turbulent flow activity. As the bottom flow changed because of the surface freshwater discharge, the fluctuation of seaward and landward flow over a short period of time caused the stress to increase in the bottom boundary layer which was linked to an increase in TKE. The Reynolds stresses showed a symmetric structure in their vertical distribution and their magnitude rapidly increased during the freshwater discharge period with a maximum magnitude of 2 <sup>×</sup> 10−<sup>4</sup> m2/s2. The TKE structure had the shape of an upside-down half-bell with an increase from the bed to 0.1–0.2 mab. Like the stresses, the magnitude of TKE was proportional to the amount of freshwater discharge and reached over <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>m</sup>2/s2 for the discharged period.

The SSC estimated from the acoustic backscatter signal of the four beams of the ADCP displayed high values in the surface and bottom layers. The high SSC in the surface happened after or during the freshwater discharge, indicating that the main source of suspended sediment was the upstream river water because there was not any other source of sediment during low tide and there was no such peak value of SSC when there was no freshwater discharged. The rapid increase of SSC in the bottom also happened after the surface freshwater discharge and it gradually propagated to the upper layer, and since it happened at the same time or after the sudden increase of the bottom TKE, the SSC near the bed was related to the intensified bottom turbulent flow that was the result of surface freshwater discharge due to the opening of the dam gate.

#### **5. Conclusions**

Using ADCP measurements in the altered Yeongsan estuary, we examined the flow in the water column and above the bed resulting from opening of the dam gate and the release of water with different properties (salinity, temperature, flow rate). The freshwater discharge was responsible for intensifying the two-layer circulation of offshore surface flow and landward bottom flow. The rapid and strong seaward surface flow affected the bottom flow structure, leading to sudden directional changes in both the horizontal and vertical flow. Responding to the freshwater discharge, the bottom turbulence also intensified rapidly, which in turn resuspended a large amount of sediment from the bed. The results of this study indicate that the surface freshwater discharge due to opening of the estuarine dam gate affects the behavior of water flow, bottom turbulence, and sediment transport in the altered Yeongsan estuary. Finally, it should be noted that since many estuarine dams have been constructed rapidly all over the world, for example, in the estuaries of the Senegal River, the Rhine-Muse rivers, and the Murray-Darling rivers, and that the construction of new estuarine dams is also under consideration, the results of this study could provide valuable insight into morphological change in estuarine environments with man-induced discrete freshwater discharges.

**Author Contributions:** Conceptualization, methodology, resources, writing—review and editing, G.-h.L.; formal analysis, investigation, validation, writing—original draft preparation, K.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Inha University Research Grant, 63142-01" and by the National Institute of Meteorological Sciences of the Korea Meteorological Administration project titled "Development of Marine Meteorology Monitoring and next-generation Ocean Forecasting System (KMA2018-00420).

**Acknowledgments:** The authors would like to thank anonymous reviewers for their constructive comments to make the high quality research.

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

#### **References**


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

## *Article* **Turbulence Characteristics before and after Scour Upstream of a Scaled-Down Bridge Pier Model**

**Seung Oh Lee <sup>1</sup> and Seung Ho Hong 2,\***


Received: 19 August 2019; Accepted: 9 September 2019; Published: 12 September 2019

**Abstract:** Bridge pier scour is one of the main causes of bridge failure and a major factor that contributes to the total construction and maintenance costs of bridge. Recently, because of unexpected high water during extreme hydrologic events, the resilience and security of hydraulic infrastructure with respect to the scour protection measure along a river reach has become a more immediate topic for river engineering society. Although numerous studies have been conducted to suggest pier scour estimation formulas, understanding of turbulence characteristics which is dominant driver of sediment transport around a pier foundation is still questionable. Thus, to understand near bed turbulence characteristics and resulting sediment transport around a pier, hydraulic laboratory experiments were conducted in a prismatic rectangular flume using scale-down bridge pier models. Three-dimensional velocities and turbulent intensities before and after scour were measured with Acoustic Doppler Velocimeter (ADV), and the results were compared/analyzed using the best available tools and current knowledge gained from recent studies. The results show that the mean flow variable is not enough to explain complex turbulent flow field around the pier leading to the maximum scour because of unsteady flows. Furthermore, results of quadrant analysis of velocity measurements just upstream of the pier in the horseshoe vortex region show significant differences before and after scour.

**Keywords:** bridge pier; horseshoe vortex; Physical hydraulic modeling; quadrant analysis; Scour and Velocity field

#### **1. Introduction**

Bridge pier scour is one type of local scour caused by sediment transport that is driven by local flow structure; therefore, it is necessary to be acquainted with the flow structure and the related scour mechanisms around the bridge pier. In general, the local flow structure around a bridge pier is composed of downflow at the upstream face of the pier in the vertical plane; the horseshoe vortex system that wraps around the base of the pier, which is the primary contributor to local scour upstream of the pier; the bow wave near the free surface on the upstream face of the pier; and the wake vortex system at the rear of the bridge pier that extends over the flow depth [1–3]. These features greatly complicate the understanding of the local flow structures [4,5], and the comprehensive effect of those complex flow structures is to increase the local sediment transport leading to additional local scour around the bridge pier [6,7].

When flow approaches a bridge pier, the velocity becomes zero on the upstream face of the bridge pier. Due to the strong adverse pressure gradient imposed by the bridge pier in the streamwise approach flow direction, the boundary layer separates upstream of the bridge pier. In the separated region, several vortices are consecutively developed, and subsequently stretched around the base of the bridge pier, giving rise to what is called a horseshoe vortex system. The primary horseshoe vortices rotate in the same sense as the approach boundary layer vorticity, but secondary vortices have the opposite rotation to preserve streamline topology. During the time that the horseshoe vortex nearest the bridge pier is decreasing in size due to stretching, a newer and younger secondary separation vortex is induced upstream of the primary vortex. The size and strength of the secondary horseshoe vortex increases with time while the size of the primary vortex continues to be reduced by stretching. At some time, the secondary and smaller vortices merge with the primary horseshoe vortex or leapfrog it to strengthen the primary horseshoe vortex, which is finally stretched completely around the bridge pier temporarily stabilizing the flow. Subsequently, instability occurs and the primary vortex forms again. The process is then repeated irregularly [8,9].

In addition to the theoretical descriptions of local vortex structures, to find more general relationships between the vortex system around a bridge pier and the resulting sediment transport, many researchers have adopted analytical, experimental, as well as computational approaches. Melville [10] has observed that the size and circulation of the horseshoe vortex increases rapidly, and the velocity near the bottom of the hole decreases as the scour hole is enlarged. According to Melville, magnitude of the downflow seems to be directly associated with the rate of scour. However, Baker [11] argued that because the size of the scour hole is much larger than the size of the vortex core on a flat-bed before scouring, such a calculation would predict wrongly that the circulation increases as scour proceeds. In his study, Baker [11] used the horseshoe vortex core circulation to derive an equation to predict the scour depth. Baker [11] assumed that the horseshoe vortex strength in the scour hole can be equal to that on a flat-bed as scour depth develops and suggested the equation with respect to pier width and free stream velocity. Later, Nakagawa and Suzuki [12] assumed that the scale and strength of the primary horseshoe vortex are constant during evolution of the local scour hole and suggested a scour prediction equation with using the stochastic nature of particle movement, originally developed by Einstein [13]. At similar time, Qadar [14] experimentally hypothesized the maximum scour depth to be a function of the initial vortex strength which is composed of vortex size and stream velocity. Based on the experimental results, an envelope curve was proposed to show the relationship between scour depth and initial vortex strength. Qadar [14] quantitatively evaluated his results by comparing with laboratory and field data.

After Qadar [14], Kothyari et al. [15] explored the diameter of primary vortex using a regression analysis of experimental data on flow around a cylindrical pier. They mentioned that the diameter of the primary vortex is dependent on the bridge opening width in comparison to the size of pier diameter. Also, Ram [16] expressed an equation for initial diameter of the horseshoe vortex and, in his findings, the initial diameter of the vortex decreases with increasing the pier Reynolds number (*Re*b) represented by a pier width as a characteristics length. However, Muzzammil and Gangadhariah [17] found that the relative vortex size, which is the ratio of vortex diameter and pier width, is weakly influenced by the pier Reynolds number for higher values on a rigid flat bed. This means that vortex size is only dependent on the pier width for higher values of the pier Reynolds number (*Re*<sup>b</sup> >104). Based on analytical models relating scour depth to horseshoe vortex properties, Muzzammil and Gangadhariah [17] proposed that equilibrium scour depth is a function of the horseshoe vortex size, tangential velocity, and vortex strength in the scour hole. They found that the mean size of the horseshoe vortex is ~20% of the cylindrical pier diameter, and the vortex tangential velocity is ~50% of the mean velocity of approach flow for 104 <sup>≤</sup> *Reb* <sup>≤</sup> 1.4 <sup>×</sup> 105 at fixed flat-bed conditions. The size of the vortex is assumed to be independent of the sediment mobility.

More recently, Dey and Raikar [7] found that the flow and turbulence intensities in the horseshoe vortex in a developing scour hole are reasonably comparable with those in before scour. Similar results can be found in other studies conducted with large-eddy simulations (LES) and detached-eddy simulation (DES) [18–20]. However, based on the experimental studies using particle image velocimetry (PIV), Guan et al. [21] argued that the size of the main vortex responsible for the maximum scour depth upstream of the pier increases with increasing scour depth.

In clear-water scour regime, as the scour depth increases, the shear stress beneath the horseshoe vortex is reduced until the shear stress becomes less than the critical shear stress and sediment movement ceases in the scour hole. As explained in the previous paragraphs, even if the role, size, and strength of turbulence leading to pier scour have been investigated qualitatively, the coherent turbulent characteristics and associated bursting events for sediment transport is still an active area of interest. So far, several of methods presented in literatures for describing the horseshoe vortex properties have not considered turbulent characteristics or unsteadiness of the horseshoe vortex; but, in fact, the bed materials around a bridge pier move irregularly in time even if the approach flow is steady. It is found in the literature that the strength and size of the horseshoe vortex are closely related to the pier geometry and the approach flow velocity upstream of the bridge while the scour depth cannot be accurately predicted unless a clear relationship between the large-scale unsteadiness of the horseshoe vortex and sediment size is presented quantitatively. One of the important objectives of this study is investigating turbulence characteristics and sequential occurrence of bursting turbulence events before and after scour in the process of particle entrainment in front of a pier. Thus, to comprehensively attack the objectives, hydraulic laboratory experiments were conducted in a flume using two different scaled-down bridge pier models. Visual observations by high speed camera as well as three-dimensional velocities and turbulent intensities before and after scour measured with ADV were analyzed using the best available tools and current knowledge.

#### **2. Methodology**

#### *2.1. Experimental Setup*

As shown in Figure 1a,b, in previous studies [22–24], laboratory experiments were conducted using various scaled hydraulic model of the Chattahoochee River bridge at Cornelia, Georgia, and Flint River bridge at Bainbridge, Georgia, USA, respectively, including the full river bathymetry. The previous experimental studies successfully explored the effect of sediment size on pier scour depth at different geometric model scales and, also based on the large number of experimental and field investigations/comparisons, improved local scour formulas were suggested. Furthermore, strategies of deciding sediment size for scaled-down hydraulic modeling was also proposed. After those model studies, each pier bent were carefully removed from the laboratory flume. Figure 1c,d show example of bridge pier bent model of Chattahoochee River bridge and Flint River bridge, respectively, with individual scales for this study. Pier bent consists of four rectangular concrete columns and rectangular concrete footings for Chattahoochee River bridge, as shown in Figure 1c. However, for Flint River bridge, two square concrete pier columns are placed on large stepped square concrete footings.

The removed pier bent model shown in Figure 1c,d was re-built inside a 1.1 m wide by 24.4 m long glass-sided tilting flume to investigate detailed turbulence and flow characteristics in front of the pier before and after scour. The approach section of the pier model was 15.0 m long followed by a working mobile bed section with a length of 3.0 m in which the pier model was placed. The length of the approach section was decided based on the findings from other researches [25–29] to ensure a fully developed approach turbulent flow and turbulent boundary layer at the bridge pier. After the working mobile bed section, additional 3 m long sediment trap section downstream of the pier model was placed. For the current set of experiments, instead of using actual cross section shape and river geometry used in the previous experiments, the channel was constructed to have a straight alignment rather than meandering, and rectangular shape of cross section was maintained through the entire flume to find more general features of flow and turbulence fields. The flume was filled with 0.53 mm bed sediment for each experiments and carefully leveled to the elevations established by the field measurements in each site. Figure 2a, b show experimental setup before running scour experiment for Chattahoochee River bridge model and Flint River bridge model, respectively.

**Figure 1.** Hydraulic laboratory model in previous laboratory studies and bridge pier model for this study: (**a**,**c**) 1:23 scaled model for Chattahoochee River bridge and (**b**,**d**) 1:33 scaled model for Flint River bridge.

**Figure 2.** Pier model in the flume: (**a**) Chattahoochee River bridge model and (**b**) Flint River bridge model.

#### *2.2. Experimental Procedure*

The experimental campaign consists of two scenarios: moveable bed experiment and fixed-bed experiment. First, moveable bed experiments were conducted to investigate variance of hydraulic parameters affecting scour depth over time. The flume was slowly filled with water to saturate the sand. After complete saturation, the required discharge (uncertainty of <sup>±</sup>2.8 <sup>×</sup> 10−<sup>4</sup> m3/s) was established with the flow depth adjusted well above the desired value. Then, the flow depth was gradually decreased by changing the height of the tailgate until the target approach flow depth was obtained. During this time, the point gage (uncertainty of ± 1 *mm*) was used to monitor the flow depth. Once the target flowrate and the flow depth had been reached, scour continued for 2 to 3 days until equilibrium was achieved. The equilibrium was defined when the increment of scour depth is less than 5% of

the bride pier diameter during 24 h. During the scouring process, instantaneous point velocities and turbulence quantities were measured by ADV in front of the pier. Furthermore, temporal change of bed elevations were measured periodically using ADV temporarily positioned for a moment above the point of scouring. At the end of scouring (equilibrium state), the velocity flow field was measured throughout the test section both in the near field next to the pier and in the far field at relative elevations for the comparison of turbulence and flow characteristics between after scour and before scour. During the velocity measurements, ADV sampling frequency was chosen to be 25 Hz with a duration of at least 2 min and perhaps as much as five minutes depending on the magnitude of turbulence at each measuring location. The correlation values in these measurements were greater than 80% and the Signal Noise Ratio (SNR) was greater than 15. The phase-space despiking algorithm was also employed to remove any spikes in the time record caused by aliasing of the Doppler signal which sometimes occurs near a boundary. More detailed filtering protocol can be found in Lee and Sturm [22] and Hong et al. [26]. After the completion of each experiment, the final bed elevations were measured using the ADV and the point gage.

When the moveable bed experiments were completed, the entire moveable bed was re-leveled and fixed by spraying polyurethane. During fixed-bed experiments, the same flow conditions as those in moveable bed experiments were reproduced, and velocities and turbulence quantities were measured in the same way as in the moveable bed experiments to investigate the effect of initial flow parameters responsible for the scour. In addition to the measurements of initial flow parameters, flow visualization experiments were also conducted. A kaolinite suspension was used as a tracer upstream of the bridge pier to show the flow structure around the bridge pier and also to detect the frequency of the horseshoe vortex system immediately upstream of the pier. The tracer was transferred from a conical tank by an electric pump operating at a maximum flowrate of 0.003 m3/s. The kaolinite suspension was mixed to achieve a concentration of 1.0 mg/cm<sup>3</sup> in the tank. The flow rate of tracer was adjusted with the aid of rotameter to produce a released velocity that was the same as the open channel mean velocity. As the tracer was released at a constant rate, a high speed video camera (30 FPS) was used to capture the unsteady dynamics of the swirl of the horseshoe vortex as it amplified and partially collapsed in size.

#### **3. Results and Discussion**

The experimental conditions have been summarized in Table 1, where *Q* is the total discharge, *b* is the width of bridge pier, *y*<sup>1</sup> is approach section water depth, *V*<sup>1</sup> is approach section velocity, *Teq* is time to the equilibrium scour, and *d*s is the equilibrium scour depth in front of the bridge pier.


**Table 1.** Summary of measured experimental conditions.

CR: Chattahoochee River bridge model and FR: Flint River bridge model.

#### *3.1. Velocity Field*

The velocity fields around the pier bent were measured for the Chattahoochee River model and Flint River model in both the fixed-bed (before scouring) and moveable-bed (after scouring in equilibrium). Figure 3a,c shows the velocity fields of "before scour" for Run 1 and Run 3, respectively. The longitudinal distance (*x*) and lateral distance (*y*) were normalized with width of corresponding bridge pier model and the near field vectors measured at 40 percent of the approach flow depth were normalized by measured approach velocity. Higher velocities were shown on both sides of the first

pier, where deeper scour occurred as the flow curved around the pier bent. In the wake zone, the mean velocities became smaller than those in the outer region. The velocity defected in this region gradually recovered in the downstream direction. The magnitude of mean velocities upstream of the first pier along the centerline became smaller approaching the pier stagnation line due to the existence of the pier. Figure 3b,d shows the combination of scour depth contours and the near-field velocities measured at 40 percent of the approach flow depth under the same flow condition as in "before scouring" but at the completion of scour for Run 1 and Run 3, respectively. The near-field velocity distributions were very close to being symmetric with respect to the centerline of the bridge pier bent as shown in Figure 3 for both runs. The characteristic decrease in magnitude was observed in near field velocity around the pier bents when comparing results before and after scour. The maximum relative difference in magnitude was approximately 30–40 percent for both cases in the vicinity of the first pier on the right-hand side. Interestingly, the maximum scour depth occurred at the nose of front pier for Run 1, as expected, however for the case of Run 3, the maximum scour occurred between the two piers with a high degree of symmetry on the left and right sides as shown in Figure 3d. For the Flint River bridge pier bent, the pier columns are placed on large stepped footings and the footings are already exposed at the beginning of scour as shown in Figure 2b. Because the footing intercepts the downflow along the nose of pier bent which feeds the size and strength of horseshoe vortex, the amount of local scour depth in front of the first pier was reduced [30–32].

**Figure 3.** Normalized mean velocity vectors: Measured at 40 percent of the approach depth for experimental Run 1 before and after scour in (**a**,**b**), respectively, and Run 3 before and after scour in (**c**,**d**), respectively.

#### *3.2. Temporal Variation of Flow and Turbulence Characteristics Upstream of the Bridge Pier*

The measured mean velocity and turbulence kinetic energy (TKE) fields have been usually used to validate three-dimensional numerical models. Even though the simulated velocity profiles at several locations were in good agreement compared with the laboratory experiments, quantitative connections with the scour depth are difficult to make. In addition, results from a three-dimensional numerical model show that the maximum mean shear stress on a fixed bed does not correspond with the maximum depth of scour hole in front of the piers [33]. Furthermore, as shown in Figure 3, the horizontal mean (time-averaged) point velocity vector plots did not indicate large changes in the velocity field with scour development when comparing the before- and after-scour conditions. It is concluded that the details of the horseshoe vortex itself must be investigated further to understand the development of the scour hole in front of the pier rather than the general-mean turbulence and flow characteristics of the near field.

Thus, in an effort to better understand the relationship between the flow field and the resulting pier scour over the development of scour hole, bed elevations, and three-dimensional velocity components as well as turbulence intensities upstream of the bridge pier were measured intermittently at two points during the scour process to capture the temporal variation of flow and turbulence characteristics as the scour hole developed. The flow was continued for the measurements so that the measurements could be completed in a time duration that was short (approximately 2–3 min) in comparison to the rate of scour hole development. The two measured points were located horizontally at a distance of one pier width upstream of the bridge pier in the streamwise direction because the size of horseshoe vortex is comparable with the size of pier width. As shown in Figure 4, the vertical location of one point (*LOC*1) was in the scour hole itself where it was varied to maintain a constant vertical displacement above the bed as the scour hole deepened with time. The other point (*LOC*2) was fixed above *LOC*1 but at a constant elevation that was close to the initial bed elevation before scour began as also shown in Figure 4. Sufficient clearance between the bridge pier and the ADV probe was required in order to place the ADV probe without bumping the pier or disturbing the horseshoe vortex system upstream of the bridge pier. Thus, the Chattahoochee River bridge model was chosen for this experiment based on the physical size and simple footing shape compared to the Flint River model.

**Figure 4.** Schematic diagram of the locations for measuring the temporal variation of flow characteristics.

The streamwise (*U*) and vertical (*W*) time-averaged velocity profiles normalized by approach mean velocity at *LOC*1 and *LOC*2 are shown in Figure 5. Both velocity profiles at *LOC*1 fluctuated slightly with time, but the values remained close to zero as the scour depth increased over time which is shown on the secondary vertical axis in Figure 5a. The averaged values became close to zero at *LOC*1 since it is located near the bed in the separation zone. However, the velocity profiles at *LOC*2 fluctuated significantly with time as shown in Figure 5b because the vertical position of *LOC*2 moved into the highly turbulent region associated with the horseshoe vortex during development of the scour hole. The fluctuations in the streamwise velocity seem to be associated with the intermittent fluctuations in the scour hole depth that result from the collapse of the sides of the hole, followed by further scouring as the hole enlarges.

The turbulence intensity in both the streamwise (*u*') and vertical (*w*') directions at *LOC*1 increased during the initial stage of scour development and then became smaller with time as shown in Figure 6a. The large fluctuations of streamwise turbulence intensity during the initial stages of scour development

is due to the unsteadiness of the location of the separation point upstream of the pier. As the scour depth changed rapidly during the initial stage, the vertical turbulence intensity became approximately four times larger than at later stages of scour hole development, indicating that the contribution of vertical turbulence intensity to the rate of scour development was very significant during the initial stage. Figure 6b shows that the turbulence intensity in both vertical and streamwise directions at *LOC*2 approached the same constant value as the scour hole developed. The vertical turbulence intensity at *LOC*2 was relatively small at the beginning of the scour process, but then increased significantly to a value greater than that at *LOC*1. This is due to the movement of the vertical location of *LOC*2 into the highly turbulent region associated with the horseshoe vortex.

**Figure 5.** Velocity profile over time at (**a**) *LOC* 1 and (**b**) *LOC* 2 for Run 1.

**Figure 6.** Temporal variations of turbulence intensity profile at (**a**) *LOC* 1 and (**b**) *LOC* 2 for Run 1.

#### *3.3. Flow Characteristics Upstream of the Bridge Pier*

While in the previous section, the long-term temporal development was considered in relationship to the variance of velocity and turbulence characteristics in the scour hole, this section investigates the short-term transient behavior of the flow immediately upstream of the bridge pier on a fixed flat bed in the region of flow separation and the horseshoe vortex. The primary horseshoe vortex is unsteady with the formation of a system of secondary vortices, and the secondary vortices are quasi-periodically combined with the primary vortex increasing its size and strength depending on the degree of stretching around the pier [8,9]. Thus, the primary horseshoe vortex oscillates in position and size in an irregular shift between two modes of behavior which are expanding and contracting over time. Those oscillation of horseshoe vortex in its size and strength in front of the pier was captured by high-speed camera in Figure 7 in this experiment, and the time difference (Δ*t*) between Figure 7a,d was related to the frequency of horseshoe vortex [22,34]. As shown in Figure 7, because of those two modes, the instantaneous velocity time series when measured near the bed close to the pier alternately exhibits periods of positive streamwise velocity towards the pier followed by negative streamwise

velocity away from the pier. The result is a bimodal velocity distribution, first described by Devenport and Simpson [34], and further studied experimentally and numerically [18,22,35]. Those two alternate states of the horseshoe vortex are bistable with the contracted mode occurring approximately 20–30% of the time. Furthermore, they found that the shape, relative size, and distance between the two peaks of the bimodal probability density function of instantaneous velocities are not permanent and stable, but instead vary with the position of the velocity measurement.

**Figure 7.** Flow visualization of horseshoe vortex with tracer injection in experimental Run 4: sequential snap shots from (**a**–**d**).

With respect to the turbulence characteristics just upstream of the pier in the horseshoe vortex region, quadrant analysis was also used to further characterize the turbulent events associated with the large-scale unsteadiness of the horseshoe vortex. Quadrant analysis was employed in this study by examining the joint probability density function of the streamwise and vertical components of fluctuating velocity, denoted as *u*' and *w*', respectively.

The four quadrants shown in Figure 8 correspond to four types of turbulent events which are defined as: I. outward interactions; II. ejections or bursts; III. inward interactions; and IV. sweeps that characterize the individual turbulent velocity measurements. The Reynolds stress is generally produced by all four types of events. The first quadrant event, I, is characterized by outward motion of high-speed fluid, with *u*' > 0 and *w*' > 0; the second quadrant event, II, is identified by outward motion of low-speed fluid, with *u*' < 0 and *w*' > 0, which is usually called ejection or bursts; the third quadrant event, III, is associated with inward motion of low-speed fluid, with *u*' < 0 and *w*' < 0 ; and finally, the fourth quadrant event, IV, represents the motion of high-speed fluid toward the bed, with *u*' > 0 and *w*' < 0, and it is called sweeps.

**Figure 8.** Schematic of the plate for quadrant analysis [36].

The relationships between the fluid motions and particle transfer near the wall can be elucidated through quadrant analysis. The ejections and sweeps contribute positively to the bed shear stress since *uw* is the flux of forward momentum to the bed, whereas the outward and inward interactions contribute negatively to the bed shear stress. The presence of a sweep corresponds to a local increase of the shear stress at the bed whereas the occurrence of an ejection corresponds to a local decrease of the shear stress at the bed. Therefore, the coherent sweep and ejection events appear to be responsible for transferring particles toward and away from the bed [37]. However, in nonuniform flows such as the wake region downstream of a backward-facing step, Keshavarzi et al. [36] have shown that turbulent events with the same shear stress contribute to different sediment transport rates due to the frequency structure of the turbulent events. While the horseshoe vortex is not expected to have the same turbulence structure as a turbulent wake, nevertheless, it does have the property of intermittency that is an important contributor to the scour process.

Among the velocity measurements made in experimental Run 2 and Run 4, two time series of velocity data at the nose of the pier were selected to conduct the quadrant analysis and to investigate the change of flow characteristics of the before-scour case and after-scour case. Figure 9 shows the results of Run 2. The time series analyzed in Figure 9a was measured for the fixed bed condition while the time series analyzed in Figure 9b was measured inside the scour hole under the same flow conditions at approximately the same distance above the bed. The joint probability density function in both figures are normalized to the peak values of 1.0, and it is necessary to multiply the contoured values by the scaling factors in the middle of the top of each figure to produce probability densities. The selected locations were both in the region of the horseshoe vortex and separated region upstream of the pier.

The results of quadrant analysis show a significant difference between the before-scour case and after-scour case in Figure 9a,b. The turbulent events for the before-scour condition were dominated by bursts and sweeps with a bimodal joint frequency distribution at the elevation of *z* = 0.17 *b*. In Figure 9a before scour, sweeps events had a higher probability of occurrence, and for both types of events, the values of *u*' were greater than those of *w*' at the maximum probability of occurrence. Run 4 shows similar probability density patterns as in Figure 9. The bursts and sweeps are the primary forcing function for creating the scour hole because they both represent positive shear stress; however, the more important characteristic is the irregular oscillation between the two types of events as the horseshoe vortex alternately expands and contracts as the separation point moves back and forth, as shown in Figure 7. The sediment particles are lifted and entrained in an intermittent fashion as explained in Lee and Sturm [22]. In the equilibrium scour hole in Figure 9b, the turbulent events no longer display the bimodal distribution with all four types of events becoming approximately equally likely. The magnitude and frequency of all events are significantly affected by the change of the flow inside the scour hole. There is no effective event of the velocity fluctuations for altering the shear stress or moving the sediment out of the equilibrium scour hole. Although the existence of the bimodal

distribution depends on the measuring location, it is significant that it disappears near the bed after scour and at the same distance above the bed as for the before-scour case.

**Figure 9.** Joint probability density function of *u*' and *w*' for (**a**) before-scour case and (**b**) after-scour case measured at *x*/*b* = −0.33 and *z*/*b* = 0.17 for Run 2.

As a further indication of the differences in the turbulence properties as the scour hole develops, the integral time scale was computed for both the before-scour and after-scour time series associated with Figure 9. The integral time scale is defined as a measure of the time over which a velocity component is dependent on its past values and a rough measure of the time interval over which a fluctuating velocity component is highly correlated with itself; it is obtained by integration of the measured autocorrelation distribution over time and a measure of the memory of the process [38]. As given in Table 2, the integral time scale in the vertical direction for the before-scour case is considerably higher than for the after-scour case; whereas, in the streamwise direction, the time scale is the same order of magnitude for both cases. The vertical fluctuation in the before-scour case is a significant contributor to the processes of suspending and eroding sediment from the bed. Accordingly, a longer integral time scale of the velocity fluctuations gives more transport than a shorter integral time scale when the rate of sediment transport increases very rapidly at the initial stage of scour development.


**Table 2.** Comparison of integral time scales for before and after scour cases.

#### **4. Summary and Conclusions**

Although local scour around bridge foundations have been extensively studied for several decades, there still remain problems because of difficulties in visualization and understanding the complex flow structure leading to scour. Thus, in this study, laboratory experiments were conducted with scaled-down bridge pier models for more complete descriptions of flow characteristics of the horseshoe vortex system around the complex bridge pier. During the experiments, velocities and turbulence intensities as well as the bed elevations before and after scour were measured by ADV. Furthermore, a simple visualization technique was used to capture the unsteadiness of horseshoe vortex. The results shows that horizontal velocity vector plot comparisons between before-scour and after-scour conditions measured at a certain relative height above the bed was not enough to explain the complex scour mechanism because the values were all temporally averaged ones, and so did not show the effects of the large-scale unsteadiness of the horseshoe vortex system upstream of a bridge pier. Further

investigation using the probability distribution of instantaneous velocity components and quadrant analysis of velocity fluctuations in the horizontal and vertical directions shows that the vertical and streamwise velocity components exhibited a bimodal probability distribution before scour near the bed upstream of the pier where horseshoe vortex is responsible for sediment transport. The streamwise and vertical velocity fluctuations were observed to be dominated by sweeps and bursts, both of which contribute positively to the bed shear stress and exhibit a bimodal joint frequency distribution. The bursts and sweeps are the primary forcing function for creating the scour hole at initial stage because they both represent positive shear stress. They are the result of irregular oscillation of the horseshoe vortex between two preferred states as the reverse flow near the bed in front of the pier either extends upstream or retreats to a point closer to the pier. As a result, the sediment particles are lifted and entrained in an intermittent fashion. After scour is complete, at equilibrium stage, the turbulent events no longer display a bimodal distribution near the bottom of the scour hole. Furthermore, a quantitative physical connection is made between the scour depth and the large-scale unsteadiness of the horseshoe vortex system in front of the pier by comparing time scales for lifting of the sediment particle, and subsequent entrainment and transport of the particle out of the scour hole.

Even if this study provides additional insights to better understand the local pier scouring process and the relationship between scour depth and the horseshoe vortex, the local flow structures in the field are affected by additional parameters, such as vertical and lateral flow contraction, unsteadiness of discharge during the passage of a flood event, and their interaction. Furthermore, in reality, scouring often happens under live-bed scour conditions with infilling of the scour hole as the flood recedes which is not well understood over a long time series because it is greatly affected by the conditions required to generally mobilize the entire bed and is further complicated by the movement of bed forms through the bridge section. Thus, additional well-designed physical model and field measurements as well as numerical simulations are required to investigate these effects, including the modeling of flow contraction, realistic hydrographs, and bed forms, such as dunes and ripples.

**Author Contributions:** S.O.L. designed/fabricated the scaled-down bridge model, and conducted the flume experiments; S.O.L. and S.H.H. provided descriptions and motivation, and analyzed the laboratory data; S.H.H. wrote the paper; all authors participated in final review and editing of the paper.

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

**Acknowledgments:** This work was partially supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MIST) (NRF-2017R1A2B2011990). Also, Seung Ho Hong was supported by West Virginia University for this study.

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

#### **References**


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

## *Article* **Turbulent Flow Field around Horizontal Cylinders with Scour Hole**

#### **Nadia Penna \*, Francesco Coscarella and Roberto Gaudio**

Dipartimento di Ingegneria Civile, Università della Calabria, Via Pietro Bucci, 87036 Arcavacata, Rende CS, Italy; francesco.coscarella@unical.it (F.C.); gaudio@unical.it (R.G.)

**\*** Correspondence: nadia.penna@unical.it; Tel.: +39-0984-496-553

Received: 21 November 2019; Accepted: 31 December 2019; Published: 2 January 2020

**Abstract:** This study presents the results of an experimental investigation on the flow-structure interactions at scoured horizontal cylinders, varying the gap between the cylinder and the bed surface. A 2D Particle Image Velocimetry (PIV) system was used to measure the flow field in a vertical plane at the end of the scouring process. Instantaneous and ensemble-averaged velocity and vorticity fields, viscous and Reynolds stresses, and ensemble-averaged turbulence indicators were calculated. Longitudinal bed profiles were measured at the equilibrium. The results revealed that suspended and laid on cylinders behave differently from half-buried cylinders if subjected to the same hydraulic conditions. In the latter case, vortex shedding downstream of the cylinder is suppressed by the presence of the bed surface that causes an asymmetry in the development of the vortices. This implies that strong turbulent mixing processes occur downstream of the uncovered cylinders, whereas in the case of half-buried cylinders they are confined within the scour hole.

**Keywords:** horizontal cylinder; turbulence structures; scour

#### **1. Introduction**

Although flow around a cylinder is one of the classical subjects of fluid dynamics, few investigations have focused on the analysis of the turbulence structures of a steady flow at scoured horizontal cylinders. This topic is relevant in the hydraulic field because it is encountered in many engineering applications, such as pipelines suspended, laid on or half-buried installed across mobile riverbeds, where a small depth-to-diameter ratio is most relevant and scouring occurs under unidirectional current [1]. In all these cases, the 3D flow field is extremely complicated due to the separation and the creation of multiple vortices. However, the complexity is further exaggerated owing to the dynamic interaction between the flow and the movable bed [2]. Erosion may occur around the pipelines, causing a higher gap between it and the bed surface and, therefore, compromising their safety. Several research findings can be found in literatures for estimating the local scour around underwater pipelines under unidirectional current [3]. Some of these studies have shown that the scour depth under unidirectional current is always higher than that under pure wave action or the combined effect of wave and current with the same bottom shear stress [3]. Accurate estimates of the scour depth are important because flow-induced oscillation by wake-vortex shedding may provoke fatigue failure of the pipeline itself [4], which is subjected to additional unsteady forces such as lift and drag.

The first investigation for estimating scour depth at submarine pipeline was conducted by Chao and Hennessy [5], who proposed an analytical method for estimating the maximum scour depth under pure current condition. The use of this method was supported by other authors e.g., [6–8]. The main drawback of the method is the use of a potential flow theory in deriving the solution. In real flow, the fluid is not inviscid, and separation occurs at the rear of the pipe. This phenomenon affects the flow conditions [3]. Later, Kjeldsen et al. [9] proposed an equation that implies that the scour

depth only depends on flow velocity and pipe diameter, but excludes the effect of flow depth and grain size. Bijker and Leeuwestein [10] stated that the scour depth depends on the undisturbed flow velocity, pipe diameter, flow depth, height of pipe above bed level, and grain size. They stated that the principal cause of erosion is a local increase in transport capacity of the water passing a pipeline, while deposition occurs where this capacity decreases [3]. Ibrahim and Nalluri [11] proposed two empirical equations relating the scour depths to flow parameters for both clear-water and live-bed conditions. The equations were derived purely from curve fitting technique. Furthermore, they stated that the grain size has no influence on the scour depth, apart from the indirect influence on the value of the critical velocity. Hansen et al. [12] obtained a relationship between the bed velocity in the scour hole and the two geometric quantities *e*/*D* and *ds*/*D*, where *e* is the gap between the underside of the cylinder and the original bed level, *D* is the cylinder diameter and *ds* the scour depth. They stated that the flume width would also influence this relationship. Mao [13] examined the scour profiles below pipelines under different flow velocities and observed that, for *ds*/*D* < 1 *ds*, is a weak function of the flow Shields parameter. He also identified two cases of the scour process: (1) jet period, which decides the maximum scour depth; (2) wake period, which decides the location of maximum scour depth. Maza [14] proposed a graphical solution for the estimation of *ds* that is a function of the initial gap-pipe diameter ratio and the flow Froude number, while Moncada and Aguirre [15] gave an empirical equation of *ds* assuming a similar functional representation used by Maza [14]. Chiew [16] identified that piping plays a dominant role in initiating scour at submarine pipelines. Later, Chiew [3] proposed an empirical function relating the flow depth ratio, *h*/*D*, with the gap-flow rate ratio, which can be used to determine the amount of gap flow through the scour hole at equilibrium conditions. Li and Cheng [17] used the finite difference method to solve the Laplace equation of velocity potential and a boundary adjustment technique to calculate the scour profiles below pipelines, while Brørs [18] used the finite element method to simulate the scour profiles below pipelines. Dey and Singh [4] examined their experimental results to describe the influence of various parameters on the equilibrium scour depth, that is flow depth, sediment gradation, different shaped cross-sections of pipes. Lately, Mohr et al. [19] related the rate of scour beneath a pipeline to the fundamental erosion properties of the sediment, namely the transport rate along the bed and the true erosion rate of the sediment. These arguments lead to two new empirical formulas that may be used to predict the time scale of the scour process beneath subsea pipelines. Note that these previous researches concentrate on soil scour around fixed pipelines. More recently, Gao et al. [20] simulated experimentally the current-induced sand scour around a vibrating pipeline to further investigate the mechanism of the coupling effects between pipe vibration and sand scour.

At the same time, several studies have been carried out both numerically and experimentally on the flow field analysis considering a flat bed. A state-of-art review of the research on the cylinder-bed surface interactions exposed to currents was conducted by Fredsøe [1]. For example, Bearman and Zdravkovich [21] carried out an experimental study on the flow around a cylinder lying horizontally at various elevation above a plane bed surface. They demonstrated that regular vortex shedding was suppressed for all gaps less than about 0.3*D*. For gaps greater than 0.3*D*, the Strouhal number, that is the ratio between inertial forces due to the unsteadiness of the flow and expresses the oscillating flow mechanisms, was found to be remarkably constant and the only influence of the plate on vortex shedding was to make it a more highly tuned process as the gap was reduced. Later Zdravkovich [22] studied in detail the flow separation from a flat plate induced by a circular cylinder. He demonstrated that when the cylinder was placed above the bed surface, the downstream separation region contained two counter-rotating vortices separated by stagnant fluid. However, no regular vortex shedding was observed in this configuration, whereas he reported that in the case of turbulent boundary layer vortex shedding occurred for a gap of 0.2*D* and for the laminar boundary layer at a gap of 0.3*D*. Lei et al. [23] analyzed the hydrodynamic forces and vortex shedding of a cylinder at different locations in the boundary layer. They proposed a quantitative method for identifying the vortex shedding suppression point. Their observations showed that the vortex shedding is suppressed at a gap of about 0.2–0.3*D*, depending on the thickness of the boundary layer. This critical gap decreases as the thickness of the boundary layer increases. Price et al. [24] performed visualization studies on the flow-cylinder interactions in the boundary layer. Specifically, they distinguished four distinct regions: (i) for gaps < 0.125*D*, the gap flow is suppressed or extremely weak, and separation of the boundary layer occurs both upstream and downstream of the cylinder. Although there is no regular vortex shedding, there is a periodicity associated with the outer shear-layer; (ii) in the region between 0.125*D* and 0.5*D*, the flow is very similar to that for very small gaps, except that there is a pronounced pairing between the inner shear-layer shed from the cylinder and the wall boundary layer; (iii) the region between 0.5*D* and 0.75*D* is characterized by the onset of vortex shedding from the cylinder; (iv) for the fourth region, with a gap greater than *D*, there is no separation of the wall boundary layer, either upstream to or downstream of the cylinder. However, downstream of the obstacle, alternate vortex shedding from the cylinder affects the wall boundary layer. Hatipoglu and Avci [25] studied experimentally and numerically the flow around a horizontal cylinder mounted and partially buried, showing that the lengths of the separation regions near the upstream and downstream of the cylinder decreased with the increasing burial ratio. More recently, Akoz et al. [26] studied quantitatively the flow characteristics of the circular cylinder laid on a fixed surface, by using the Particle Image Velocimetry (PIV) technique. The main purpose was to reveal the mechanisms of vortical flow structures that are mostly responsible for scour and burial processes. They demonstrated that the intersection of the bed surface and cylinder enhances the burial mechanisms hydrodynamically even in wake flow regions. Furthermore, it was shown that the wake flow region is shortened in size in the longitudinal direction as a function of the Reynolds number. More recently, Arslan et al. [27] studied for four different submergence ratios the 3D unsteady flow around a rectangular cylinder with the large eddy simulation (LES) turbulence model. To conduct the analysis, they used the experimental data obtained by Malavasi and Guadagnini [28]. As a result, they characterized the behavior of vortex structures generated by separated flow and the hydrodynamic forces acting on the semi-submerged structure.

Few are the studies in which the turbulence flow structures are analyzed focusing in the neighborhood of the cylinder considering a mobile bed. The first studies on this topic were those of Kjeldsen et al. [9] and Mao [13], who examined the flow around a pipe placed over a scour hole, Jensen et al. [29] investigated experimentally the flow around a pipeline placed initially on a flat bed at five characteristic stages of the scour process in currents. They found that the mean flow field and turbulence around and the force on a pipeline undergo considerable changes, as the scour below the pipeline develops in time and space. Moreover, vortex shedding comes into existence at early stages of scour process, first in a somewhat premature form caused by the close proximity of the dune formed downstream the pipe as a result of deposition of scoured material. As the dune moves away from the pipe, the vortex shedding gradually reaches a stage which resembles the shedding process of a free cylinder. Starting from this experimental campaign, Smith and Foster [30] examined the flow physics around the pipeline (considering also the bed scouring phenomenon) using numerical modelling. They compared the mean horizontal and vertical velocities and the wake characteristics measured by Jensen et al. [29] to those obtained by applying a computational fluid dynamics (CFD) model (Flow3D) using both a two-equation *k*-ε model and a Smagorinsky LES turbulence closure scheme for five stages of the scour process.

Thus, the main purpose of the study was to identify the flow structures in shallow water condition developed once the bed scour reached the equilibrium, varying the gap between the cylinder and the bed surface. Specifically, three different conditions were reproduced in a laboratory flume: a suspended cylinder, a laid on cylinder and a half-buried cylinder. The analysis was focused on the velocity measurements in correspondence of the horizontal cylinder and on the determination of the instantaneous and ensemble-averaged velocity fields, the instantaneous and ensemble-averaged vorticity fields, viscous and Reynolds stresses, and the ensemble-averaged turbulence indicators. In addition, the equilibrium longitudinal bed profiles were discussed in order to better comprehend the effect induced by the flow-structure interactions.

#### **2. Experimental Set-Up and Procedure**

The experimental campaign was performed at the Laboratorio "Grandi Modelli Idraulici" (GMI, Università della Calabria, Italy). A recirculating flow channel (9.6 m long, 0.485 m wide, 0.5 m deep) was used in this study. The flume side walls were made of glass in order to visualize the flow. The inlet of the flume comprised of a stilling tank, an uphill slipway, and honeycombs (having a diameter of 10 mm) to dampen the flow disturbances at the entry. The test section was located at 7.33 m downstream of the flume entrance and it was 0.165 m long. The flow depth was regulated by a downstream tailgate. To collect the outflow, a tank equipped with a calibrated Thomson weir to measure the flow discharge was attached downstream of the tailgate.

The bed was constituted by very coarse sand having a median sediment size *d*<sup>50</sup> of 1.53 mm (0.06 mm < *d* < 2 mm, where *d* is the size of sediments) and geometric standard deviation σ*<sup>g</sup>* [=(*d*84/*d*16) 0.5] of 1.24, where *d*<sup>16</sup> and *d*<sup>84</sup> are the 16% and 84% (by weight) finer sizes of sediments, respectively. The sediment density was ρ*<sup>s</sup>* = 2680 kg/m3. To prepare the bed, sediments were initially spread within the flume and screeded to create a bed with a longitudinal slope *S*<sup>0</sup> of 0.1%.

A horizontal cylinder of 30 mm in diameter made of Plexiglas was set up at the center of the test section. The cylinder was fixed to the flume walls and kept normal to the flow direction. Three different gaps, respectively *e*1, *e*<sup>2</sup> and *e*3, between the cylinder and the bed surface were considered, as represented in Figure 1 and reported in Table 1: (i) suspended cylinder; (ii) cylinder laid on the bed; (iii) partially buried cylinder. To ensure direct comparisons of the results, the bed surface was leveled at the start of each experiment.

**Figure 1.** Schematic representation of the geometrical conditions of the experimental runs, where *ei* (*i* = 1, 2, 3) is the gap between the lower edge of the cylinder and the bed surface.


**Table 1.** Geometric and hydraulic parameters of the experimental runs.

The runs were carried out under the same hydraulic conditions. Specifically, the initial flow depth *h* was about 0.074 m, as measured immediately upstream to the cylinder from the sediment crest. Therefore, the aspect ratio *B*/*h* was greater than 5 (*B* being the flume width) and the effect of the sidewalls is negligible [31] in our experimental set-up. However, the ratio *h*/*D* was kept constant and equal to 2.5. Thus, being *h*/*D* < 5 [4], the shallowness effect is present in all the runs and, therefore, the resulting scouring process is due to the cylinder-bed surface interactions at different gaps in shallow water condition. The threshold flow velocity, *Uc*, was determined using Neill's empirical formula [32] and was equal to 0.37 m/s, higher than the average flow velocity *U* = *Q*/(*Bh*) ≈ 0.28 m/s, *Q* being the flow discharge. This indicates a clear-water condition. The flow Reynolds number *Re* (=4*Uh*/υ, where <sup>υ</sup> <sup>=</sup> 1.04 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>2/s at 18.4 ◦C is the kinematic water viscosity), the flow Froude number *Fr* [=*U*/(*gh*) 0.5, where *g* is the gravitational acceleration], the cylinder Reynolds number *Rec* (=*UD*/υ, where *D* is the cylinder diameter) and the other hydraulic parameters are listed in Table 1 for each experimental run. Specifically, upstream of the cylinder (where the flow was not affected), the shear velocity *u*\* was estimated by extending the Reynolds shear stress profile linearly to the maximum gravel crest as −*uw*- 0.5 *z*=*zc* (where *u* and *w* are the temporal velocity fluctuations in the streamwise and vertical directions, respectively, and *zc* is the maximum gravel crest elevation). Hence, the Reynolds number of the sediments *Re*\* was calculated as *u*\*ε/υ, where ε is the equivalent sand roughness height ≈ 2*d*<sup>50</sup> [33].

During all the experimental runs bed scouring occurred. It was found that the equilibrium scour depth was reached within 24 h for all the three experimental runs. This is in line with the findings of Chiew [3]. Dey and Singh [4] demonstrated that the time to reach the equilibrium scour below a pipeline is shorter than that at an abutment, since it is due to strong pressurized flow and not to primary or horseshoe vortex. Also, Mao [13] observed equilibrium scour below a cylinder in just over 3 h.

A 2D PIV system manufactured by TSI was used to measure the flow field. It consisted of a 12 bit CCD camera (Nikon, city, country, 50 mm F1.8 lens, 2048 <sup>×</sup> 2048 pixel2 resolution, 15.2 <sup>×</sup> 15.2 mm2 sensor size) and a double pulse Nd:YAG laser with a frame rate of 15 Hz and a pulse energy of 200 mJ at a wavelength of 532 mm. Titanium dioxide powder, having a diameter of 3 μm and a mass density of 4.26 kg/m3, was used as a tracer during flow measurements. The test section was illuminated with the laser to capture the movement of the tracer particles through the camera placed parallel to the laser sheet. The flow measurements were taken in correspondence of the horizontal cylinder. In order to measure the flow field, 3000 pairs of images were captured over a period of about 200 s. The inter-frame time between two laser pulses was set equal to 1400 μs, which means that, considering *U* ≈ 0.28 m/s, we were able to measure only eddies greater than 0.39 mm. The actual size of the eddy is imposed by the spatial resolution of the PIV measurements [34]. The field of view was 165 <sup>×</sup> 165 mm2 with an interrogation area of 32 <sup>×</sup> 32 pixel2. A 50% overlap of the interrogation areas was employed to increase the spatial resolution of the measurement [35] to about 1.3 mm (which corresponds to 12.4 pixel/mm). The Insight 4G-2DTR software was used during the acquisition phase and to process the resulting data.

The precision in the flow velocity measurements with a PIV system depends to a great extent on the errors introduced by the sub-pixel estimator in the cross-correlation. This is known as "peak-locking" and is a bias error in the PIV that occurs when the particle images are too small, causing the particle image displacement to be biased toward integer value. The error is estimated to be 10% of the particle image diameter, which is the diameter in pixels of the particle as seen through the camera [36]. The mean particle image diameter in the present case is about 1 pixel, and a typical displacement between the cross-correlation image pairs is 5 pixels in the main flow direction. Therefore, the estimated random error of the measured velocity vector in each interrogation area is about 2% for the streamwise velocity component.

However, as *a priori* method, in order to reduce peak-locking, through a trial and error procedure we set an aperture number *f* # equal to 11, in order to ensure that all the tracer particles in the light sheet (having a thickness of 2 mm) were in focus. Nevertheless, following the procedure described in Padhi et al. [34], we used also and *a posteriori* method pre-processing all the images by using a filter (already available in the Insight 4G-2DTR software), which optimized the particle image diameter with respect to the peak estimator, in order to improve the signal-to-noise ratio.

The flow was stopped at the end of the measurement phase and, subsequently, the flume was slowly emptied without any disturbance for the bed scour by using a bottom outlet located 1 m downstream of the test section. For each run, the bed surface was acquired with the photogrammetry technique. Photogrammetry is used in different fields such as topographic mapping, architecture, engineering, cultural heritage and geology e.g., [37–39]. Following the procedure described in Penna et al. [40], a Nikon D5200 camera was used, equipped with a Nikkor 18–55 mm f/3.5–5.6 G VR lens. The resulting 3D point cloud was at first transformed into an unstructured triangular mesh by using the software PhotoScan (Agisoft, St. Petersburg, Russia). Then a structured grid was extracted using the commercial software Rhinoceros (McNeel & Associates, Seattle, WA, USA), with a spatial resolution δ*l* = 5 mm in both the streamwise and spanwise directions.

#### **3. Results and Discussion**

#### *3.1. Longitudinal Bed Profiles*

The centerline longitudinal bed profiles at equilibrium, for each experimental run, are shown in Figure 2. Here, the origin of the *x*-axis (streamwise direction) is set in correspondence of the cylinder center, whereas the *z*-axis indicates the vertical direction. Specifically, the horizontal and vertical axes were made dimensionless by dividing them by the cylinder diameter. Furthermore, *z*/*D* = 0 is the bed surface at the start of each experiment. These longitudinal profiles were extracted from the bed surface models derived from the photogrammetry technique.

**Figure 2.** Centerline longitudinal bed profiles at equilibrium.

The data clearly highlight how the gap between the cylinder and the bed surface influences the scour formation and, thus, the bed profile. A steady flow that impacts on a cylinder laid on the bed causes a deeper scour hole than in the case in which the cylinder is suspended. Specifically, it was 1.5*D* = 4.6 cm and 1.9*D* = 5.7 cm for Run 1 and Run 2, respectively. This means that the gap in Run 1 allows the flow to pass with less disturbance to the sand bed. It is also evident that the maximum equilibrium scour depth *ds* occurs closer to the cylinder in the second case than in Run 1. The distance between the center of the cylinder and the maximum equilibrium scour depth, *xs*, is 5*D* = 15 cm for Run 1 and 3.5*D* = 10.5 cm for Run 2. The behavior can be explained by analyzing the vortex system induced by the fluid-structure interaction, which is described in the following sections. Furthermore, it is possible to note that the resulting scour hole involves the area upstream to and beneath the cylinder, owing to a particular process also known as tunnel erosion [41], that is ascribed to seepage processes [2,42] and occurs when the ratio between the flow depth and the cylinder diameter is less than 3.5 [16]. The vortices due to the obstacle led to instabilities to the particles, moving them away. The consequent erosion process is due to the increased velocity underneath the cylinder.

A different behavior can be noted considering Run 3 in which the cylinder is half-buried: tunnel erosion does not occur and the scour hole develops downstream of the cylinder itself. Hence, in this case, the shape of the bed profile is comparable to that formed in the presence of a bed sill subjected to a steady flow. This will be further explained analyzing the flow field around the cylinder in the following sections. Interestingly, even if *ds* is the smallest value (*ds* = 1.5*D* = 4.4 cm) if compared to the other two cases, the scour length *ls* is about 0.3*D* = 0.9 m and, therefore, the scouring process affects a larger amount of sediments than in Run 1 and Run 2, for which it is almost 0.22*D* = 0.65 m. Note that the small hump visible in Figure 2 at *x*/*D* = 0 is the profile of the cylinder itself.

The obtained bed profile at equilibrium for *e* = 0 (Run 2) was compared with the experimental data of Chiew [3], Gao et al. [20] and Dey and Singh [4]. To better understand the different behavior that can be observed in Figure 3, Table 2 shows the experimental conditions adopted for each run. The bed profiles of Gao et al. [20] and Dey and Singh [4] were determined in deep-water condition (*h*/*D* > 5), whereas those related to the studies conducted by Chiew [3] and of the present work were in shallow-water condition (*h*/*D* < 5). Therefore, it is easy to note that for shallow-water condition a deeper and longer scour hole is formed, leaving aside for the moment considerations about the other parameters involved. Furthermore, in the case of deep-water condition the scour hole develops to great extent below the cylinder, whereas for shallow-water the maximum scour depth is reached downstream of it. By comparing the bed profile of Run 2 and that of Gao et al. [20], which are characterized by similar cylinder diameter and flow velocity, the higher water depth of Gao et al. [20] produced a maximum scour depth of about 1.8 cm, whereas that obtained in the present study was about 5.7 cm, although *d*<sup>50</sup> was 4 times less than 1.53 mm of Run 2. However, it must be pointed out that also the considered equilibrium time (*te*) can influence the observed discrepancies between the bed profiles. Based on results obtained from Chiew [3], only 50–70% of the equilibrium scour depth is reached in three to four hours of testing. Finally, by comparing the bed profile of Run 2 and that of Dey and Singh [4], which are characterized by a similar sediment size, it is possible to note that the maximum scour depth is of the same order, despite the fact that *D* and *U* are higher than those of the present study. This can again be ascribed to both the water depth and also to the equilibrium time.

**Figure 3.** Comparison between centerline longitudinal scour profiles around a horizontal cylinder with *e* = 0.

**Table 2.** Geometric and hydraulic parameters of the literature experimental runs for *e* = 0.


#### *3.2. Instantaneous and Ensemble-Averaged Velocity Fields*

A sequence of four consecutive instantaneous flow fields at different times *t* for all the experimental runs are shown in Figure 4. The velocity has a magnitude <sup>|</sup>**u**<sup>|</sup> <sup>=</sup> *u*<sup>2</sup> + *w*<sup>2</sup> 0.5, where *<sup>u</sup>* and *<sup>w</sup>* are the instantaneous streamwise and vertical velocity components, respectively. The velocity was made dimensionless by dividing it by *U*. The origin of the *x*-axis refers to the start of the test section acquired with the PIV system. Both the axes were made dimensionless by dividing by *D*. Furthermore, Figure 4 shows the streamlines for each velocity field. Note that invalid vectors were removed from the instantaneous velocity fields (blank areas in the colormaps), as well as those below and along the edge of the cylinder since they corresponded to an area poorly illuminated by the laser. Vectors identified as spurious were replaced by the ensemble-averaged values when computing the average over time.

**Figure 4.** Dimensionless instantaneous velocity fields at different times *t* on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

From the analysis of Figure 4, it is evident that upstream to the cylinder, at the beginning of the study area, the flow velocities, in all the three runs, increase with the vertical distance *z*, as usually observed in an open-channel flow, and the streamlines are nearly horizontal. Thus, in both Run 1 and Run 2, the incoming flow approaches the upstream face of the cylinder and causes flow separation. The main flow is divided into two flows: one part is oriented towards the bed surface (this is the main cause of scouring) and the other part is oriented towards the free surface [26], giving birth to two separation zones. Since 300 < *Rec* < 3 <sup>×</sup> 105 (Table 1), the boundary layer over the cylinder surface is laminar. For Run 2 the separation zone at the free-stream side of the cylinder moves upstream with respect to that of Run 1. This is in agreement with previous findings [26,29,43]. A different behavior occurs for Run 3: the obstruction represented by the cylinder creates an adverse pressure gradient resulting in the separation of flow lines forming a vortex on the upstream front of the cylinder. Separation of flow also takes place at the free-stream side of the cylinder. Reverting to Run 1 and Run 2, a pair of vortices in the wake of the cylinder is noticeable, causing vortex shedding. The vortices are shed alternately at both sides of the cylinder at a certain frequency. Specifically, following the classification of Sumer and Fredsøe [43], the generated wake is completely turbulent, whereas the vortex-shedding frequency, *f*, can be derived from the Strouhal number:

$$St = \frac{f \times D}{U} \tag{1}$$

The Strouhal number in proximity of a wall may change with respect to the case in which no wall is present. Shedding frequency tends to increase (yet slightly) as the gap ratio decreases [44]. However, the flow field in this study was analyzed at the equilibrium phase, that is when the scouring phenomena was already exhausted. This means that the gap between the cylinder and the bed is 4.27 cm and 3.43 cm for Run 1 and Run 2, respectively (Figure 2). Therefore, since the gap ratio *ei*/*D* is higher than 1 for Run 1 and Run 2 (1.42 and 1.14, respectively), it is possible to assume that in these cases the Strouhal number in proximity of the wall is equal to the Strouhal number for wall-free cylinder [44,45]. Note that the vortex shedding occurs in a region where the scour hole is deeper than that at the cylinder itself (*e*1/*D* > 1.42 and *e*2/*D* > 1.14 for Run 1 and Run 2, respectively). Therefore, the wall represented by the mobile bed has no effects on the Strouhal number, even if concave and not flat. Considering that for a smooth circular cylinder *St* remains practically constant at the value of 0.2 for 300 < *Rec* < 3 <sup>×</sup> 105 [46], the vortex-shedding frequency can be determined as *St* <sup>×</sup> *U*/*D*, thus 1.87 s−<sup>1</sup> for Run 1 and Run 2. This means that vortex shedding occurs about twice per second. Figure 4 helps in visualizing the formation and destruction of the vortex system, from *t* = 0.13 s to 0.53 s, with *t* = 0 at the beginning of the experimental measurements. The vortex originated from the free-stream side of the cylinder (e.g., Figure 4b at *t* = 0.13 s) is strong enough to draw the opposing vortex (at the cylinder wall-side) downstream across the wake region (e.g., Figure 4b at *t* = 0.27 s). Note that the streamlines coming from above are oriented in the clockwise direction, while the others in the anti-clockwise direction. The opposite sign will then cut off further supply of vorticity to the upper vortex from its boundary layer [43], and therefore this vortex is shed and the transported downstream by the flow (e.g., Figure 4b at *t* = 0.40 s). Thus, a new vortex forms at the free-stream side of the cylinder. However, the lower vortex is now stronger than this new upper vortex and, therefore, this will lead to its shedding (e.g., Figure 4b at *t* = 0.53 s). Instead, for Run 3 a reverse roller is formed in the wake region. Vortex shedding is suppressed by the presence of the bed surface that causes an asymmetry in the development of the vortices [43].

To analyze the ensemble-averaged flow field in a spatial flow domain, Figure 5 shows the streamlines and the contours of the velocity having a magnitude **u** <sup>=</sup> *u*<sup>2</sup> + *w*<sup>2</sup> 0.5, where *<sup>u</sup>* and *<sup>w</sup>* are the ensemble-averaged streamwise and vertical velocity components, respectively. The ensemble-averaged velocity was made dimensionless by dividing it by *U*.

**Figure 5.** Dimensionless ensemble-averaged velocity fields on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

The analysis of Figure 5 recalls the separation of the main flow upstream to the cylinder in both Run 1 and Run 2, described on the basis of the instantaneous velocity fields of Figure 4. It is apparent that the increase of the gap between the cylinder and the bed surface (Run 1) causes the formation of two opposite vortices in the wake of the cylinder shorter and more compressed than those of Run 2. This behavior may be attributed to the fact that in Run 1 the cylinder is close to the water surface, pushing down the upper vortex. In the meantime, a greater amount of flow passes below the cylinder and, therefore, tends to expand immediately downstream of it, compressing the vortex system in the wake region. The velocity at the underside of the cylinder is slightly lower than the one at the top side. This is in agreement with the observations of Jensen et al. [29]. The ensemble-averaged velocity field of Run 3, instead, shows the roller vortices upstream to and downstream of the cylinder in the scour hole, as discussed before. This led to a flow acceleration along the water depth above the recirculation zone.

Figure 6 shows the vertical profiles of the ensemble-averaged streamwise flow velocity at different distances upstream to and downstream of the cylinder. Again, it is illustrated how the velocity profiles are not influenced by the horizontal cylinder at the beginning of the area of interest, in all the three experimental runs. For Run 3 only, negative values of *u* are shown upstream to the cylinder, owing to the presence of the vortex originated from the interaction between the current and the obstacle. Immediately downstream of the cylinder in Run 1 and Run 2, the streamwise flow velocity is negative, indicating a reverse flow. Increasing the distance from the cylinder, it loses intensity and the vertical profiles start recovering the undisturbed upstream condition. However, the velocity profiles are different if we refer to Run 3, in which the shear layer separates from the cylinder, creating a recirculation zone.

**Figure 6.** Vertical profiles of the ensemble-averaged streamwise flow velocity at different streamwise distances in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

#### *3.3. Instantaneous and Ensemble-Averaged Vorticity Fields*

The contours of dimensionless instantaneous vorticity ω*yD*/*U* are shown in Figure 7 at different times *t.* Here, ω*<sup>y</sup>* is the instantaneous vorticity given by ∂*u*/∂*z* − ∂*w*/∂*x*, whose positive and negative values specify clockwise and counterclockwise fluid motions, respectively. Specifically, the counterclockwise fluid motion causes the flow to accelerate, resulting in a downward transport of momentum in the downstream direction, and the clockwise fluid motion causes the flow to decelerate, resulting in an upward transport of momentum in the upstream direction [47].

**Figure 7.** Dimensionless instantaneous vorticity fields at different times *t* on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

From the analysis of Figure 7, the patterns of vorticity show that in the wake region of Run 1 and Run 2 a cluster of vorticity takes place along the shear layers. This is due to the separation induced by a pressure increase along the cylinder surface. The vortex that originates from the separation zone at the free-stream side of the cylinder has a positive sign indicating a clockwise rotation, whereas the vortex from the wall-side has a negative sign with a counterclockwise rotation. These motions are responsible for the production of the Reynolds shear stresses [48]. It is also evident that for Run 2 the steep slope of the scour hole behind the cylinder forces the shear layer originating from the lower edge of the cylinder to bend upwards, thus causing the associated lower vortex to interact with the upper one prematurely, as it was observed by Jensen et al. [29]. As regards to Run 3, the vortex downstream of the cylinder rotates with clockwise direction in the scour hole, creating a recirculation zone with decelerated fluid motion. However, in all the three cases it is found that the vortices do not change their rotational sense with time; therefore, they are recursive with time.

Figure 8 shows the streamlines and the contours of the dimensionless ensemble-averaged vorticity ω*yD*/*U*, to substantiate the effects illustrated by the instantaneous vorticity patterns. Here, ω*<sup>y</sup>* is the ensemble-averaged vorticity given by ∂*u*/∂*z* − ∂*w*/∂*x*. At the upstream end of the investigated area, the vorticity is null along the water depth in all the three experimental runs. However, in the near-bed flow clockwise vortices occur owing to the bed roughness. For Run 1 and Run 2, downstream of the cylinder, in the wake region, vorticity assumes positive values in the upper zone and negative values in the lower part, that is with clockwise and counterclockwise rotations, respectively. Instead, Run 3 shows the clockwise vortex in the scour hole, as previously discussed.

**Figure 8.** Dimensionless ensemble-averaged vorticity fields on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

#### *3.4. Viscous and Reynolds Stresses*

The variations of the dimensionless viscous stresses τ*<sup>v</sup>* in the flow domain, calculated as (ν*du*/*dz*)/*u*<sup>2</sup> <sup>∗</sup> are depicted in Figure 9. Instead, the dimensionless Reynolds shear stresses τ*uw* = −*uw*-/*u*<sup>2</sup> <sup>∗</sup> are shown in Figure 10.

By comparing Figures 9 and 10, it is evident that upstream to the cylinder the viscous stresses are negligible along the water depth and increase in the near-bed flow zone in all the experimental runs. At the same time, τ*uw* tends to assume a peak value at a certain distance from the bed surface and then starts decreasing. A different behavior is shown in Run 3, since a small region of negative values of τ*uw* is evident close to the bed surface, owing to the small vortex on the upstream front of the cylinder. Considering Run 1 and Run 2, downstream of the cylinder two regions characterized by the maximum positive and negative values of the Reynolds shear stress occur. This suggests strong turbulent mixing process in these areas, indicating a downward transport of momentum in the downstream direction and an upward transport of momentum in the upstream direction, respectively. Close to the cylinder surface, τ*<sup>v</sup>* assumes the highest positive and negative values where the shear stress is almost null.

As regards Run 3, both τ*<sup>v</sup>* and τ*uw* assume their maximum values in the scour hole, which suggests a strong turbulent mixing process that entrains high-momentum fluid into the recirculation region.

**Figure 9.** Dimensionless viscous stresses on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

**Figure 10.** Dimensionless Reynolds shear stresses on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

Finally, the dimensionless streamwise and vertical Reynolds normal stresses, expressed as σ*uu* = −*uu*-/*u*<sup>2</sup> <sup>∗</sup> and σ*ww* = −*ww*-/*u*<sup>2</sup> <sup>∗</sup>, respectively, are shown in Figures 11 and 12.

**Figure 11.** Dimensionless streamwise Reynolds normal stresses on a vertical plane in (**a**) Run 1, (**b**) Run 2, (**c**) Run 3.

**Figure 12.** Dimensionless vertical Reynolds normal stresses on a vertical plane in (**a**) Run 1, (**b**) Run 2, (**c**) Run 3.

In all the experimental runs, upstream to the cylinder, it can be seen that σ*uu* and σ*ww* assume their maximum negative values in the near-bed flow region, and then they increase with *z*. This may be attributed to the fluid mixing that occurs in the presence of bed roughness. At the same time, it can be noted that downstream of the cylinder the intense fluid mixing induced by the presence of the obstacle enhances both the fluctuations *u* and *w*- . The comparison between Run 1 and Run 2 shows that a smaller gap between the cylinder and the bed surface induces the flow area to enlarge with increased magnitude of both *u* and *w*- .

On the other hand, Run 3 shows downstream of the cylinder a decrease of σ*uu* and σ*ww* in the recirculation zone, indicating, also in this case, fluid mixing with a maximum magnitude of 0.08 m/s and 0.05 m/s for *u* and *w*- , respectively, that are less than those in Run 1 and Run 2, for which *u* and *w*are about 0.14 m/s and 0.09 m/s. This leads to a decrease of about 65% for both σ*uu* and σ*ww*.

#### *3.5. Ensemble-Averaged Turbulence Indicators*

In order to measure the local turbulence level in a flow, the use of the ensemble-averaged turbulence indicator I may represent a useful tool [47]. This parameter can be estimated as follows:

$$I = \left(\frac{\frac{2}{3}k}{lL^2}\right)^{0.5},\tag{2}$$

where *k* is the turbulent kinetic energy (TKE), expressed as:

$$k = \frac{1}{2} (\overline{\overline{u'u'}} + \overline{\overline{v'v'}} + \overline{\overline{w'w'}}) , \tag{3}$$

*v* being the fluctuations of the spanwise velocity component with respect to its ensemble-averaged value. Since the PIV system allows the measurement of only two velocity components, the term *vv*- can be approximated as 0.5 *uu*- + *ww*- [49,50]. Therefore, the TKE is expressed as:

$$k = 0.75(\overline{u'u'} + \overline{w'w'}).\tag{4}$$

The turbulence level is high if *I* > 0.5, moderate if 0.1 < *I* < 0.5 and low if *I* < 0.1 [51]. Thus, the contours of the ensemble-averaged turbulence indicator for the three experimental runs are shown in Figure 13. It is clear that for Runs 1 and 2 the magnitudes of *I* are moderate immediately downstream of the cylinder (*I* ≈ 0.3), where the mean velocity is very low and a reverse flow takes place. In the wake region, the fluctuations in the streamwise and vertical directions cause the enhancement of the turbulence level, but it can be still classified as moderate (*I* ≈ 0.5 for Run 1 and *I* ≈ 0.4 for Run 2). However, the area characterized by the enhancement of the turbulence level is greater in Run 2 than

in Run 1, since in this case it interests the entire flow depth. This can be ascribed to the fact that the increase of the gap between the cylinder and the bed surface (Run 1) causes the formation of two opposite vortices shorter and more compressed than those of Run 2. This also leads to higher *I* value in Run 1 than in Run 2. As regards Run 3, the turbulence level reaches a maximum value of about 0.4 in the scour hole downstream of the cylinder and, therefore, it can be classified as a moderate one. Nevertheless, in the rest of the flow domain, the turbulence level is basically low and it decreases as *z* increases, implying that the roughness effects and the influence of the cylinder diminishes owing to the damping in *u* and *w*- , which is visible also in Figures 11 and 12.

**Figure 13.** Ensemble-averaged turbulence level on a vertical plane in (**a**) Run 1, (**b**) Run 2 and (**c**) Run 3.

#### **4. Conclusions**

This study examines the turbulence characteristics at scoured horizontal cylinders in different condition of submergence, in order to deepen the knowledge investigating the flow-structure interactions in the presence of a mobile bed in shallow-water condition. Three experimental runs were performed considering a suspended cylinder, a laid on cylinder and a half-buried cylinder. Specifically, bed scouring occurred in all the runs and the resulting bed configurations were analyzed in order to better comprehend the flow dynamics.

For the first two conditions, at the end of the scouring phase, it was demonstrated that the incoming flow approaches the upstream face of the cylinder and causes flow separation. Thus, the main flow is divided into two flows: one part is oriented towards the bed surface, whereas the other part is oriented towards the free surface. The velocity at the underside of the cylinder is slightly lower than the one at the top side, as observed by Jensen et al. [29] in their experiments. This latter is also responsible for the formation of the scour hole that involves the area upstream to and beneath the cylinder owing to tunnel erosion. This phenomenon is more pronounced when the cylinder is laid on the bed. In the wake zone, for the suspended cylinder a greater amount of flow passes below it. Therefore, this flow tends to expand immediately downstream of the cylinder, limiting the extension of the induced vortex system. The analysis of the vorticity field revealed that in the wake region vorticity assumes positive values in the upper zone and negative values in the lower part, that is with clockwise and counterclockwise rotations, respectively. This is in accordance with the results obtained by Bearman and Zdravkovich [21], Price et al. [24] and others, since it was demonstrated that vortex shedding occurs at gaps higher than 0.3*D*. The strong turbulent mixing process in these areas is highlighted by the Reynolds shear stress distribution, indicating a downward transport of momentum in the downstream direction and an upward transport of momentum in the upstream direction, respectively. To support these observations, the turbulence level was calculated and it was found that fluctuations in the streamwise and vertical directions cause the enhancement of the turbulence level immediately downstream of the obstacle, reaching a value within the range 0.3–0.5, indicating a moderate turbulence level.

A different behavior was observed when the cylinder was initially half-buried. Tunnel erosion did not occur and, as a consequence, the obstruction represented by the cylinder created an adverse pressure gradient resulting in the separation of flow lines forming a vortex on the upstream front of the cylinder. Separation of flow also took place at the free-stream side of the cylinder. A reverse roller was formed in the wake region. This was responsible for the local scour downstream of the cylinder. The vorticity pattern showed the presence of the clockwise vortex in the scour hole with a strong mixing process that entrains high-momentum fluid into the recirculation region. It was demonstrated that here the turbulence level can be classified as moderate (*I* ≈ 0.4).

Additional work in the future is required to take into account other variables that may affect the investigated turbulence structures (e.g., cylinder diameter, shape of the pipe, flow discharge, water depth), their development during the scouring process and also other in-depth statistical analyses (such as the third-order statistics to provide a very accurate measure of the TKE dissipation rate [52] or correlation functions to analyze the turbulent coherent structures [53]).

This study represents an advancement in the current understanding of the flow-structure interaction at scoured horizontal cylinders subjected to currents. It was demonstrated that bed changes should be considered in the analysis of the flow field and in the prediction of the forces acting on the cylinder rather than considering a plane boundary, because the effects are very different and depend also on the gap between the cylinder and the bed itself.

**Author Contributions:** Conceptualization, N.P. and R.G.; Data curation, N.P. and F.C.; Formal analysis, N.P. and F.C.; Methodology, N.P., F.C. and R.G.; Writing—Original Draft Preparation, N.P.; Writing—Review and Editing, N.P., F.C. and R.G.; Supervision, R.G. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The authors would like to thank Antonio Leuzzi for his valuable work during the performance of the experimental Runs.

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

#### **References**


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