**1. Introduction**

Ice jams in rivers are a common ice phenomenon in cold regions in winter that affect both hydraulic conditions and boundary conditions of watercourses. Sometimes, under the same flow condition, the presence of an ice cover in a river leads to a higher water level, which can induce ice flood disasters. The presence of bridge piers in a river narrows the flow cross section, changes the flow state of water around the piers, and forms a horseshoe vortex, which often results in a local scour around bridge piers. As reported by researchers, the existence of bridge piers in rivers increases the possibility of the formation of ice jams [1,2], and the formation and development of an ice jam aggravate the local scour of bridge piers [3]. Therefore, the mutual influence of flow conditions, local scour, and ice jam evolution makes relevant problems more complex.

Recently, some research has been conducted to study the hydraulic interaction between bridge piers and ice jams. The results of laboratory experiments showed that the existence of bridge piers results in an increase in the ice transport capacity near bridge piers and a

**Citation:** Hu, H.; Wang, J.; Cheng, T.; Hou, Z.; Sui, J. Channel Bed Deformation and Ice Jam Evolution around Bridge Piers. *Water* **2022**, *14*, 1766. https://doi.org/10.3390/ w14111766

Academic Editor: Giuseppe Oliveto

Received: 19 May 2022 Accepted: 27 May 2022 Published: 31 May 2022

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

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

decrease in the ice jam thickness [2]. Wang et al. [4,5] carried out an experimental study on ice jam blockage at bridge piers, which showed that ice blockage was affected by the ice cube size, spanning distance between piers, and flow conditions. Based on experiments conducted in a curved flume in a laboratory, it was found that the bridge pier located at the apex of the river bend had less influence on ice jam formation than the one installed in a straight section of a river (the junction of two bends) [6]. The results of laboratory experiments regarding the variation of water level around piers with different shapes under an ice jam showed that the increment of water level in the present of a cylindrical pier was the smallest, while it was the largest for a rectangular pier, and the regression formula for calculating the increase in water level caused by an equilibrated ice jam was given [7]. It was found that the ice accumulation process around bridge piers depends on the pier size, ice discharge rate, flow discharge, and Froude number. For a flow with a smaller Froude number and a larger ice discharge rate, with the presence of a smaller pier in the channel, the evolution process of an ice jam is similar to that without the presence of a pier [8]. There are very few studies on numerical simulation. Cheng et al. [9] used a discrete element model to simulate the ice floe accumulating process around bridge piers. In addition, Yang et al. [10] used a material point method; Istrati et al. [11] used a finite element method; Salciarini et al. [12] used a discrete element method, and Hasanpour et al. [13] used a smoothed particle hydrodynamics method to study tsunami-induced debris impacts on bridge structures, respectively, which has a certain reference value for the numerical simulation of hydraulic action between bridge piers and ice jams. Parola et al. [14] conducted experimental studies to develop a method for predicting debris forces on piers and superstructures. Malavasi et al. [15] studied the hydrodynamic interactions between an open flow and a bridge deck. The results of this work will provide a theoretical reference for the study of the interaction between ice floes and bridge decks under extreme hydraulic conditions.

Some researchers have done research work regarding the local scour process around circular piers under both open flow and ice-covered (instead of ice-jammed) flow conditions. In an open flow, to improve the applicability of the equation developed by Melville and Coleman in laboratory and field data, Pandey et al. [16] proposed new K-factors for a nonuniform gravel bed. Memar et al. [17] studied the influence of the size and position of circular collars from the sediment bed on serial bridge piers. Zhao et al. [18] derived one formula for calculating scour depth around a skew bridge in a curved road by solving the deflection angle of water flow and the width of the water crossing section. The existence of an ice cover changes the flow structure [19–22] and affects the development of a local scour. Compared to the adjacent position of the collar on the bed, placing the collars below the bed would reduce the scouring speed around a pier. Bacuta and Dargahi [23], Ackermann and Shen [24], and Hains and Zabilansky [25,26] studied the local scour at cylindrical bridge piers under ice-covered flow conditions and found that compared to open flow conditions, the local scour depth is greater under an ice-covered flow condition. Munteanu [27] performed experiments under four different boundary conditions (open flow, totally ice-covered, two sides partially ice-covered and one side partially ice-covered) and pointed out that the local scour process around the bridge pier was the most intense under two sides in a partial ice-covered condition, and the maximum scour depth was about 55% higher than that under an open flow condition. Christopher et al. [28] investigated the scour problem around bridge piers under different ice cover thicknesses and found that the scour depth increased with the increase in ice cover thickness. Wu et al. [29–31] investigated the local scour process under ice-covered flow conditions based on laboratory experiments and reported that the maximum scour depth increased with the increase in cover roughness. They also established empirical formulas for calculating the scour depth under an ice-covered flow condition. Wang et al. [32] analyzed the relationship between the flow intensity and dimensionless scour depth, as well as the difference in the rate of change of scour depth under an ice-covered flow condition compared to that under an open flow condition and developed an empirical equation for determining the change of

scour depth with the time under an ice-covered condition. Namaee and Sui [33–37] studied the influence of ice covers with different roughness on the local scour process around side-by-side piers. They found that compared with the smooth-covered flow condition, the near-bed flow velocity was higher, and the maximum flow velocity position was closer to the channel bed under the rough cover. The maximum scour depth around piers increased with the decrease in particle size and the increase in cover roughness and densimetric Froude number. The smaller the pier size and the larger the pier spacing distance, the weaker the horseshoe vortex around the pier and the shallower the scour hole. The grain size of the armour layer obviously affects the depth of the scour hole.

The existing research studies either focus on the influence of bridge piers on the stability and evolution of ice jams or the local scour process around piers under an ice-covered (instead of ice-jammed) flow condition. However, no research work has been reported regarding the local scour process around bridge piers during the evolution process of an ice jam considering the thickness of an ice jam and the time for achieving the equilibrium condition. In the present experiment study, the influence of a local scour at bridge piers on the evolution of ice jams was studied. The development of scour holes under different evolution modes of the ice jam was investigated. The results of the present study can provide a basis for bridge design and safety protection in practical engineering.

### **2. Materials and Methods**

### *2.1. Experimental Setup*

Laboratory experiments were carried out in a flume. The plan layout of the flume is shown in Figure 1. The length of the flume is 26.68 m with a width of 0.4 m and a depth of 1.3 m. Along this flume from upstream to downstream, 22 cross sections (CS) for measurement with an equal spacing distance of 1.2 m were set up. Between the upstream CS-4 and CS-5, an ice hopper was setup for discharging model ice particles into the flowing water. A Styrofoam panel that was 0.6 m long and 0.4 m wide was placed on the water surface between CS-20 and CS-21 for assisting with the initiation of an initial ice jam in front of the foam panel in the flume, corresponding to the phenomenon of an initial ice jam head formation in natural rivers.

Polyethylene particles were used to model ice particles. The mass density of the polyethylene particles was 0.918 g/cm<sup>3</sup> , which is nearly equal to the mass density of ice of 0.917 g/cm<sup>3</sup> in rivers. Polyethylene particles has a flat ellipsoid shape with the longest diameter of 3.5 mm (as shown in Figure 1c). Between CS-2 and CS-22, a sand bed was prepared to model a sand bed in rivers. The initial thickness of the sand bed in the flume was 10 cm. The median particle size (*d*50) of the three sands used in the experiments was 0.49 mm, 0.713 mm, and 1.04 mm, respectively. The inhomogeneity coefficient of these three sands was 2.0, 1.61, and 2.68, respectively. The mass density was 1.423 g/cm<sup>3</sup> . Three different shapes of model piers were used, namely cylindrical, rectangle, and round end-shaped piers, as shown in Figure 1. The model pier was installed in the flume center at CS-16. Cross section 3 was the control section for the initial flow depth (*H*0) and initial velocity (*V*0) and the ice discharge rate (*Q*<sup>i</sup> ) from the ice hopper. In engineering design, the ratio of the pier width to a single span between adjacent piers is 0.04~0.16 [38]. Considering the flume width of 40 cm, the pier width (*D*) was 2 cm. In addition, this research work is a conceptual study, which is not targeted at a specific engineering prototype, but mainly investigates the changes of local scours and median particle sizes of sediment on the speed of ice jam evolution and scour hole development.

**Figure 1.** (**a**,**b**) Flume layout for experiment; (**c**) polyethylene particles; (**d**) cross section of model piers. **Figure 1.** (**a**,**b**) Flume layout for experiment; (**c**) polyethylene particles; (**d**) cross section of model piers.

#### *2.2. Experimental Process 2.2. Experimental Process*


ice jam was the average of the largest and smallest values along this cross section; when an ice wave phenomenon appeared at this cross section, the thickness of the ice jam at this cross section was the average of the thickness of the wave crest and wave trough). The water level was measured at the same frequency as that for the ice jam thickness, and both were performed simultaneously until the ice jam and local scour achieved an equilibrium state. Experiments showed that the time for the scour hole to achieve an equilibrium state under an ice-jammed flow condition was related to the hydraulics and ice flow conditions and ranged from 8 h to 11 h. However, all experiments lasted 24 h to ensure that the local scour under an ice jam condition achieved an equilibrium state. the average of the largest and smallest values along this cross section; when an ice wave phenomenon appeared at this cross section, the thickness of the ice jam at this cross section was the average of the thickness of the wave crest and wave trough). The water level was measured at the same frequency as that for the ice jam thickness, and both were performed simultaneously until the ice jam and local scour achieved an equilibrium state. Experiments showed that the time for the scour hole to achieve an equilibrium state under an ice-jammed flow condition was related to the hydraulics and ice flow conditions and ranged from 8 h to 11 h. However, all experiments lasted 24 h to ensure that the local scour under an ice jam condition achieved an equilibrium state.

(6) Subsequently, we stopped adding ice particles in the flume, raised the tailgate, and gradually reduced the flow discharge to increase the water level so that ice particles under an ice jam in the flume stopped moving. To obtain a precious bathymetry of the scour hole, laboratory measurements were carefully carried out in order to not disturb the scour hole and ice jam due to the addition of probes. A point gage and a spot measuring instrument with an accuracy of 0.1 mm were used for measuring the bathymetry of the deformation of the sand bed around the pier. For each experimental run, data were collected at 26 measurement points around the bridge pier. Some variables used in this study were defined and are shown in Figure 2. (6) Subsequently, we stopped adding ice particles in the flume, raised the tailgate, and gradually reduced the flow discharge to increase the water level so that ice particles under an ice jam in the flume stopped moving. To obtain a precious bathymetry of the scour hole, laboratory measurements were carefully carried out in order to not disturb the scour hole and ice jam due to the addition of probes. A point gage and a spot measuring instrument with an accuracy of 0.1 mm were used for measuring the bathymetry of the deformation of the sand bed around the pier. For each experimental run, data were collected at 26 measurement points around the bridge pier. Some variables used in this study were defined and are shown in Figure 2.

**Figure 2.** Scour hole at the pier under an ice jam and associated variables. ((**a**): scour hole with point measurement; (**b**,**c**): points for measuring ice jam thickness). **Figure 2.** Scour hole at the pier under an ice jam and associated variables. ((**a**): scour hole with point measurement; (**b**,**c**): points for measuring ice jam thickness).

Each experimental run was considered to achieve an equilibrium state if the following conditions were met: (1) Water levels at all cross section did not change; (2) The thickness of the ice jam at any cross section did not change, namely, the amount of ice particles discharged from the ice hopper located upstream equalled that of the output at the downstream end of the ice jam; (3) The scour hole also reached an equilibrium condition, namely, not only the depth and shape of the scour hole around the pier did not change, but also the height and shape of the deposition dune downstream of the scour hole did not change. Whenever an experimental run satisfied these conditions, this experimental run was considered to reach the equilibrium state. Each experimental run was considered to achieve an equilibrium state if the following conditions were met: (1) Water levels at all cross section did not change; (2) The thickness of the ice jam at any cross section did not change, namely, the amount of ice particles discharged from the ice hopper located upstream equalled that of the output at the downstream end of the ice jam; (3) The scour hole also reached an equilibrium condition, namely, not only the depth and shape of the scour hole around the pier did not change, but also the height and shape of the deposition dune downstream of the scour hole did not change. Whenever an experimental run satisfied these conditions, this experimental run was considered to reach the equilibrium state.

The evolution of an ice jam is divided into three stages: (1) The formation of the initial ice jam phase: the ice jam develops from the downstream Styrofoam panel to upstream with the continuous constant incoming ice particles from the ice hopper (as shown in Figure 3a). (2) The ice jam thickening stage: with the continuous incoming ice particles from the ice hopper, the ice jam gains thickness from upstream to downstream since the incoming ice particles are entrained by the flowing current and submerged at the head of an ice jam and then gradually delivered downstream. This process continues before the ice jam reaches an equilibrium state (as shown in Figure 3b). (3) The ice jam equilibrium stage: the amount of ice particles coming from upstream is equal to the output amount of ice particles (the amount of ice particles collected from the pond located at the downstream of the flume) from the toe of the ice jam body; however, the appearance of The evolution of an ice jam is divided into three stages: (1) The formation of the initial ice jam phase: the ice jam develops from the downstream Styrofoam panel to upstream with the continuous constant incoming ice particles from the ice hopper (as shown in Figure 3a). (2) The ice jam thickening stage: with the continuous incoming ice particles from the ice hopper, the ice jam gains thickness from upstream to downstream since the incoming ice particles are entrained by the flowing current and submerged at the head of an ice jam and then gradually delivered downstream. This process continues before the ice jam reaches an equilibrium state (as shown in Figure 3b). (3) The ice jam equilibrium stage: the amount of ice particles coming from upstream is equal to the output amount of ice particles (the amount of ice particles collected from the pond located at the downstream of the flume) from the toe of the ice jam body; however, the appearance of waves of ice accumulation migrate from upstream to downstream, and as a consequence, the ice jam becomes thicker

whenever the crest of a wave of ice accumulation passes and thinner whenever the trough of such as wave passes (as shown in Figure 3c). In the present study, to assess the local scour around a pier, various experiments with local scour (as shown in Figure 4) and without local scour (namely, the fixed-bed experiments) were conducted, as summarized in Table 1. the ice jam becomes thicker whenever the crest of a wave of ice accumulation passes and thinner whenever the trough of such as wave passes (as shown in Figure 3c). In the present study, to assess the local scour around a pier, various experiments with local scour (as shown in Figure 4) and without local scour (namely, the fixed-bed experiments) were conducted, as summarized in Table 1. thinner whenever the trough of such as wave passes (as shown in Figure 3c). In the present study, to assess the local scour around a pier, various experiments with local scour (as shown in Figure 4) and without local scour (namely, the fixed-bed experiments) were conducted, as summarized in Table 1.

waves of ice accumulation migrate from upstream to downstream, and as a consequence, the ice jam becomes thicker whenever the crest of a wave of ice accumulation passes and

waves of ice accumulation migrate from upstream to downstream, and as a consequence,

*Water* **2022**, *14*, x FOR PEER REVIEW 6 of 20

*Water* **2022**, *14*, x FOR PEER REVIEW 6 of 20

(a) Initial ice jam phase (b) Ice jam thickening stage

(c) Ice jam equilibrium stage

**Figure 3.** Development of an ice jam at different stages.

**Figure 4.** Local scour hole around a pier. **Figure 4.** Local scour hole around a pier. **Figure 4.** Local scour hole around a pier.


**Table 1.** Summary of all experimental runs.

Pier shape factor *K<sup>ζ</sup>* refers to a table for calculating the pier shape factor and pier width in the appendix of the Code for Hydrological Specifications for Survey and Design of Highway Engineering (JTG C30-2015). The pier shape factors of different pier types involved in this study are shown in Table 2.

**Table 2.** The pier shape factors corresponding to different pier shape.


### **3. Interaction of Local Scour and Ice Jam Evolution**

*3.1. Impact of Local Scour on Ice Jam Evolution*

Figure 5 shows the vertical distribution of flow velocity at the cross section where the bridge pier is located, with and without a scour in the vicinity of the pier. One can see in Figure 5 that the difference in velocity distribution before and after an initial ice jam reaches the cross section (CS-16) where the pier is located (as shown in Figure 6b). It should be noted here that when an initial ice jam does not reach CS-16, it means that the head of an initial ice jam is located between CS-17 and CS-16 or downstream of CS-16 (as shown in Figure 6a), while the ice jam that passes CS-16 corresponds to the fact that the head of an initial ice jam is located upstream of CS-15 (as shown in Figure 6c).

(a) The ice jam did not reach the bridge pier (a) The ice jam did not reach the bridge pier

(b) The ice jam reaches the bridge pier (b) The ice jam reaches the bridge pier

(c) The ice jam passes through reach the bridge pier (c) The ice jam passes through reach the bridge pier

**Figure 6.** Pier scour for 3 cases: before the head of an initial ice jam reaching the pier (**a**); the head of an initial ice jam reaching the pier (**b**); after the head of an initial ice jam passing the pier (**c**). **Figure 6.** Pier scour for 3 cases: before the head of an initial ice jam reaching the pier (**a**); the head of an initial ice jam reaching the pier (**b**); after the head of an initial ice jam passing the pier (**c**). **Figure 6.** Pier scour for 3 cases: before the head of an initial ice jam reaching the pier (**a**); the head of an initial ice jam reaching the pier (**b**); after the head of an initial ice jam passing the pier (**c**).

As shown in Figure 5, before the head of an initial ice jam reaches and passes the pier, the flow velocity near the water surface or at the bottom of the initial ice jam at CS-16 is smaller. Also, the flow velocity near the channel bed is relative higher with the presence of a pier scour, compared to that without a local scour (fixed-bed condition), which makes the initial ice jam develop faster around the pier. When the head of an initial ice jam passes the pier (CS-16), the development of the initial ice jam is mainly influenced by the flow condition (velocity and depth) and the scour conditions (with or without a scour). As seen from Figure 7, in the presence of a pier scour, water levels upstream of the pier (CS-16) are greater after the initial ice jam passes CS-16; thus, the initial ice jam develops more rapidly toward the upstream section of the channel. The comparison of the time required to complete the formation of an initial ice jam along the entire channel with and without a local scour is shown in Table 3. Compared with the results without a local scour process around the pier, the time required to complete the formation of an initial ice jam along the channel is shorter with the presence of a local scour around the pier. Because of the local scour process around the pier, the flow velocity near the pier is smaller, which leads to a more rapid development of the initial ice jam toward upstream. Meanwhile, with the propagation of the head of an initial ice jam toward upstream, the water levels upstream of the initial ice jam head are higher, and thus the flow velocity is lower, and it takes a shorter time to complete the formation of an initial ice jam along the channel.

As shown in Figure 5, before the head of an initial ice jam reaches and passes the pier, the flow velocity near the water surface or at the bottom of the initial ice jam at CS-16 is smaller. Also, the flow velocity near the channel bed is relative higher with the presence of a pier scour, compared to that without a local scour (fixed-bed condition), which makes the initial ice jam develop faster around the pier. When the head of an initial ice jam passes the pier (CS-16), the development of the initial ice jam is mainly influenced by the flow condition (velocity and depth) and the scour conditions (with or without a scour). As seen from Figure 7, in the presence of a pier scour, water levels upstream of the pier (CS-16) are greater after the initial ice jam passes CS-16; thus, the initial ice jam develops more

*Water* **2022**, *14*, x FOR PEER REVIEW 9 of 20

rapidly toward the upstream section of the channel.

**Figure 7.** Comparison of water levels in the upstream cross section when the initial ice jam passes the pier (CS-16) with a local scour and those without a local scour. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s; **c**: *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**d**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s). **Figure 7.** Comparison of water levels in the upstream cross section when the initial ice jam passes the pier (CS-16) with a local scour and those without a local scour. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s; (**c**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**d**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s).

The comparison of the time required to complete the formation of an initial ice jam along the entire channel with and without a local scour is shown in Table 3. Compared with the results without a local scour process around the pier, the time required to complete the formation of an initial ice jam along the channel is shorter with the presence of a local scour around the pier. Because of the local scour process around the pier, the flow velocity near the pier is smaller, which leads to a more rapid development of the initial ice jam toward upstream. Meanwhile, with the propagation of the head of an initial ice jam toward upstream, the water levels upstream of the initial ice jam head are higher, and thus the flow velocity is lower, and it takes a shorter time to complete the formation of an initial ice jam along the channel.


*Water* **2022**, *14*, x FOR PEER REVIEW 10 of 20

**Table 3.** The time required to complete the development of an initial ice jam phase with and without local scour conditions. local scour conditions.

**Table 3.** The time required to complete the development of an initial ice jam phase with and without

The flow depths under an equilibrium ice jam with a local scour were compared to those without a local scour, as presented in Figure 8. At each designated cross section, under the same hydraulic condition and ice discharge rate, the flow depth under an equilibrium ice jam without a local scour process ("B") is smaller than that with a local scour process ("A"), and thus, the average flow velocity is larger, and the ice transport capacity is stronger. This means that the existence of a local scour should consume part of the energy of the flow under an ice jam, which is not conducive to the transportation of ice particles to the downstream section. The flow depths under an equilibrium ice jam with a local scour were compared to those without a local scour, as presented in Figure 8. At each designated cross section, under the same hydraulic condition and ice discharge rate, the flow depth under an equilibrium ice jam without a local scour process ("B") is smaller than that with a local scour process ("A"), and thus, the average flow velocity is larger, and the ice transport capacity is stronger. This means that the existence of a local scour should consume part of the energy of the flow under an ice jam, which is not conducive to the transportation of ice particles to the downstream section.

**Figure 8.** Comparison of water depth under an equilibrium ice jam with a local scour around the pier to those without a local scour. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s and *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s and *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s). **Figure 8.** Comparison of water depth under an equilibrium ice jam with a local scour around the pier to those without a local scour. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s and *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s and *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s).

Figure 9 shows the comparison of water levels during the state of an equilibrium ice jam with a local scour to those without a local scour around a pier. It can be seen from Figure 9 that for an equilibrium ice jam, the water level at each cross section with a local scour is greater than that without a local scour around the pier. At the same time, it can be seen that in general, with the presence of a local scour at the pier, the increase in the water level along the channel reaches between CS-2 and CS-4 is the highest. This phenomenon results from the thickening stage of an ice jam, namely, the thicker the ice jam, the higher the upstream water level and the larger the hydraulic slope, as reported by Sui et al. [39]. The numbers in Figure 9 show that when the water level for an equilibrium ice jam without the presence of a local scour is used as the benchmark value, one can clearly see the rate of increase in the water level for an equilibrium ice jam with the presence of a local scour condition. The rate of increase in water level at the cross section where the pier is located is also larger because the presence of a local scour and the pier resulted in an increase in water level compared with that at the other cross sections. On the other hand, one can see from Figure 10 that the thickness of an equilibrium ice jam at each cross section with a local scour is greater than that without a Figure 9 shows the comparison of water levels during the state of an equilibrium ice jam with a local scour to those without a local scour around a pier. It can be seen from Figure 9 that for an equilibrium ice jam, the water level at each cross section with a local scour is greater than that without a local scour around the pier. At the same time, it can be seen that in general, with the presence of a local scour at the pier, the increase in the water level along the channel reaches between CS-2 and CS-4 is the highest. This phenomenon results from the thickening stage of an ice jam, namely, the thicker the ice jam, the higher the upstream water level and the larger the hydraulic slope, as reported by Sui et al. [39]. The numbers in Figure 9 show that when the water level for an equilibrium ice jam without the presence of a local scour is used as the benchmark value, one can clearly see the rate of increase in the water level for an equilibrium ice jam with the presence of a local scour condition. The rate of increase in water level at the cross section where the pier is located is also larger because the presence of a local scour and the pier resulted in an increase in water level compared with that at the other cross sections. On the other hand, one can see from Figure 10 that the thickness of an equilibrium ice jam at each cross section with a local scour is greater than that without a local scour around a pier. Because of the existence of a local scour, the water flow will lose part of the energy through the bridge pier section, which

reduces the ice transport capacity of the water flow near the bridge pier section, and during the evolution of the ice jam, the ice particles near the pier are transported downstream less, and the thickness of the cross-section ice jam increases. While the thickness of the ice jam near the pier thickens, the upstream water level increases, thereby reducing the water flow ice transport capacity of the upstream section, and there is more ice particle accumulation in the upstream section.

Table 4 shows the comparison of the time required for an ice jam to reach an equilibrium condition with a local scour process around the pier to that without a local scour. As shown in Table 4, it takes a longer time for the ice jam to reach an equilibrium jam stage with the presence of a local scour than that without a local scour at the pier. As mentioned above, the presence of a local scour consumes the flow energy near the pier. Because of the presence of a local scour at the pier, the water depth of the flow cross section increased. During the thickening stage of an ice jam, the ice transport capacity around the pier is weakened. Thus, there are more ice particles accumulated at the bottom of the ice jam near the pier; as a consequence, it causes an increase in the thickness of the ice jam and a reduction in the water depth under the ice jam. This will affect the transport capacity of ice particles. This local scour process and ice jam evolution interact with each other, and thus, it takes a longer time for an ice jam development to achieve an equilibrated stage with the presence of a local scour at the pier. *Water* **2022**, *14*, x FOR PEER REVIEW 11 of 20 local scour around a pier. Because of the existence of a local scour, the water flow will lose part of the energy through the bridge pier section, which reduces the ice transport capacity of the water flow near the bridge pier section, and during the evolution of the ice jam, the ice particles near the pier are transported downstream less, and the thickness of the crosssection ice jam increases. While the thickness of the ice jam near the pier thickens, the upstream water level increases, thereby reducing the water flow ice transport capacity of the upstream section, and there is more ice particle accumulation in the upstream section.

**Figure 9.** Comparison of water levels during an equilibrium ice jam with a local scour around the pier to those without a local scour. The numbers in Figure 9 show that when the water level for an equilibrium ice jam without the presence of a local scour is used as the benchmark value, the rate of increase in water level for an equilibrium ice jam with the presence of a local scour condition is in %. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s; (**c**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**d**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s). **Figure 9.** Comparison of water levels during an equilibrium ice jam with a local scour around the pier to those without a local scour. The numbers in Figure 9 show that when the water level for an equilibrium ice jam without the presence of a local scour is used as the benchmark value, the rate of increase in water level for an equilibrium ice jam with the presence of a local scour condition is in %. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s; (**c**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**d**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s).

**Figure 10.** Comparison of the thickness of an equilibrium ice jam with the presence of a local scour at the pier to that without a local scour. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s; (**c**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**d**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s). **Figure 10.** Comparison of the thickness of an equilibrium ice jam with the presence of a local scour at the pier to that without a local scour. ((**a**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.17 m/s; (**b**): *H*<sup>0</sup> = 0.25 m, *V*<sup>0</sup> = 0.15 m/s; (**c**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.17 m/s; (**d**): *H*<sup>0</sup> = 0.20 m, *V*<sup>0</sup> = 0.15 m/s).


Table 4 shows the comparison of the time required for an ice jam to reach an equilibrium condition with a local scour process around the pier to that without a local scour. As shown in Table 4, it takes a longer time for the ice jam to reach an equilibrium **Table 4.** Comparison of the time required for an ice jam to reach an equilibrium condition with and without a local scour around the pier.

**Table 4.** Comparison of the time required for an ice jam to reach an equilibrium condition with and without a local scour around the pier. **Bridge Pier Shape** *V***<sup>0</sup> (m/s)** *H***<sup>0</sup> (m)** *<sup>Q</sup>***<sup>i</sup> (L/s) Time Required to Complete the Initial Ice Jam Phase (min) With Local Scour Without Local Scour With Local Scour Without Local Scour** cylindrical 0.17 0.25 0.0205 0.0205 8.5 7.0 cylindrical 0.15 0.25 0.0205 0.0205 9.0 8.0 cylindrical 0.17 0.20 0.0205 0.0205 8.0 7.0 cylindrical 0.15 0.20 0.0205 0.0205 8.5 7.5 round end-shaped 0.17 0.20 0.0205 0.0224 7.5 5.5 round end-shaped 0.16 0.20 0.0205 0.0212 8.5 6.5 As seen in Figure 11, for the same ice discharge rate (*Q*<sup>i</sup> ), the ratio of jam thickness (*Ti* ) to flow depth under an ice jam (*h*) *Ti*/*h* decreases with the increase in the flow Froude number *F<sup>r</sup>* regardless of whether a local scour appears around the pier. An increase in the *F<sup>r</sup>* means that the ice transport capacity of the flow under the ice jam increases, and more ice particles are transported downstream during the evolution of the ice jam (as shown in Figure 12). Thus, the thickness of the ice jam becomes relatively small (as shown in Figure 13). It was also found that with the same *F<sup>r</sup>* , the *Ti*/*h* is larger when there is a local scour around the pier. Because of the presence of a local scour, the flow velocity at the cross section where the pier is located decreases. Therefore, the transport capacity of the ice decreases; the ice jam evolution process slows down, and fewer ice particles are delivered to the downstream of the pier. Thus, the thickness of the ice jam increases, resulting in a greater value of *Ti*/*h* in the equilibrium state. During the thickening process of an ice jam near the pier, the upstream water level increases. Thus, the ice transfer capacity in the

upstream section is also reduced. More ice particles accumulate in the upstream section. The flow velocity under the ice jam increases and causes the local scour upstream of the pier. Thus, the *Ti*/*h* value in the upstream cross section is greater (as shown in Figure 14). In Figure 12, *F<sup>r</sup>* is the flow Froude number of the initial flow at the cross section where the pier is located. upstream section is also reduced. More ice particles accumulate in the upstream section. The flow velocity under the ice jam increases and causes the local scour upstream of the pier. Thus, the *Ti*/*h* value in the upstream cross section is greater (as shown in Figure 14). In Figure 12, *Fr* is the flow Froude number of the initial flow at the cross section where the pier is located. upstream section is also reduced. More ice particles accumulate in the upstream section. The flow velocity under the ice jam increases and causes the local scour upstream of the pier. Thus, the *Ti*/*h* value in the upstream cross section is greater (as shown in Figure 14). In Figure 12, *Fr* is the flow Froude number of the initial flow at the cross section where the pier is located.

As seen in Figure 11, for the same ice discharge rate (*Q*i), the ratio of jam thickness (*Ti*) to flow depth under an ice jam (*h*) *Ti*/*h* decreases with the increase in the flow Froude number *Fr* regardless of whether a local scour appears around the pier. An increase in the *Fr* means that the ice transport capacity of the flow under the ice jam increases, and more ice particles are transported downstream during the evolution of the ice jam (as shown in Figure 12). Thus, the thickness of the ice jam becomes relatively small (as shown in Figure 13). It was also found that with the same *Fr*, the *Ti*/*h* is larger when there is a local scour around the pier. Because of the presence of a local scour, the flow velocity at the cross section where the pier is located decreases. Therefore, the transport capacity of the ice decreases; the ice jam evolution process slows down, and fewer ice particles are delivered to the downstream of the pier. Thus, the thickness of the ice jam increases, resulting in a greater value of *Ti*/*h* in the equilibrium state. During the thickening process of an ice jam near the pier, the upstream water level increases. Thus, the ice transfer capacity in the

As seen in Figure 11, for the same ice discharge rate (*Q*i), the ratio of jam thickness (*Ti*) to flow depth under an ice jam (*h*) *Ti*/*h* decreases with the increase in the flow Froude number *Fr* regardless of whether a local scour appears around the pier. An increase in the *Fr* means that the ice transport capacity of the flow under the ice jam increases, and more ice particles are transported downstream during the evolution of the ice jam (as shown in Figure 12). Thus, the thickness of the ice jam becomes relatively small (as shown in Figure 13). It was also found that with the same *Fr*, the *Ti*/*h* is larger when there is a local scour around the pier. Because of the presence of a local scour, the flow velocity at the cross section where the pier is located decreases. Therefore, the transport capacity of the ice decreases; the ice jam evolution process slows down, and fewer ice particles are delivered to the downstream of the pier. Thus, the thickness of the ice jam increases, resulting in a greater value of *Ti*/*h* in the equilibrium state. During the thickening process of an ice jam near the pier, the upstream water level increases. Thus, the ice transfer capacity in the

*Water* **2022**, *14*, x FOR PEER REVIEW 13 of 20

*Water* **2022**, *14*, x FOR PEER REVIEW 13 of 20

**Figure 11.** Relationship between *Ti*/*h* and *Fr* with and without a local scour. **Figure 11.** Relationship between *Ti*/*h* and *Fr* with and without a local scour. **Figure 11.** Relationship between *Ti*/*h* and *Fr* with and without a local scour.

**Figure 12.** Transportation of ice particles. **Figure 12.** Transportation of ice particles. **Figure 12.** Transportation of ice particles.

**Figure 13.** Relationship between *Ti*/*h* and *Fr* with the presence of a local scour. **Figure 13.** Relationship between *Ti*/*h* and *Fr* with the presence of a local scour.

*3.2. Effect of Ice Jam Evolution on Local Scour*

the ice jam:

**Figure 14.** The ratio *Ti*/*h* at different cross sections with and without a local scour with the same *Fr.*

1. The bed shear stress changes dramatically after an ice cover forms on the water surface compared to that of open flow conditions. Thus, in the presence of an ice cover, the turbulence intensities in a scour hole increase, and the sediment transport will be affected significantly. As pointed out by Wu et al. [29–31], the rougher the ice cover, the closer the maximum velocity to the sand bed and thus, the greater the Reynolds shear stress inside the scour hole. As a consequence, the scour hole at a pier becomes deeper. When the development of an ice jam is dominated by a mechanical thickening process (an ice cover is unable to withstand the external forces acting on it and is broken up and extruded to form an ice jam) [40], the development process of a scour hole is similar to that under an open flow condition. In such a case, as the ice jam develops upstream, the thickness of the ice jam increases. When the initial ice jam reaches the pier, the thickness of the initial ice jam grows rapidly, and the flow depth decreases. Thus, the velocity under the initial ice jam at this section increases and results in an intensified scouring process around the pier with a rapid increase in the scour depth. The scouring process at the pier gradually decreases since the increase in the jam thickness gradually slows down and gradually approaches the equilibrium ice jam phase, as shown in Figure 15 (In the figure, *L*<sup>i</sup> is the length of a

scour hole, and *h*<sup>i</sup> is the depth of a scour hole at the pier, respectively).

2. When the formation of an initial ice jam is dominated by a hydraulic thickening process (ice floes and frazil ice particles are entrained by flowing water and submerged at the front of the ice cover and accumulate under the ice cover) [40], with

The results show that compared to the formation of a scour hole under an open flow condition, the formation of a scour hole under an ice-jammed flow condition can be classified into the following two different categories based on the thickening process of

**Figure 13.** Relationship between *Ti*/*h* and *Fr* with the presence of a local scour.

**Figure 14.** The ratio *Ti*/*h* at different cross sections with and without a local scour with the same *Fr.* **Figure 14.** The ratio *Ti*/*h* at different cross sections with and without a local scour with the same *Fr*.

#### *3.2. Effect of Ice Jam Evolution on Local Scour 3.2. Effect of Ice Jam Evolution on Local Scour*

The results show that compared to the formation of a scour hole under an open flow condition, the formation of a scour hole under an ice-jammed flow condition can be classified into the following two different categories based on the thickening process of the ice jam: The results show that compared to the formation of a scour hole under an open flow condition, the formation of a scour hole under an ice-jammed flow condition can be classified into the following two different categories based on the thickening process of the ice jam:


*Water* **2022**, *14*, x FOR PEER REVIEW 15 of 20

**Figure 15.** Comparison of the development of a scour hole around a pier under open flow, ice cover and ice jam conditions in Case 1. **Figure 15.** Comparison of the development of a scour hole around a pier under open flow, ice cover and ice jam conditions in Case 1. **Figure 15.** Comparison of the development of a scour hole around a pier under open flow, ice cover and ice jam conditions in Case 1.

the scouring process around the pier slows down, as shown in Figure 16.

the scouring process around the pier slows down, as shown in Figure 16.

fewer ice particles entrained and accumulating under an ice cover/initial ice jam, the thickness of the initial ice jam is much smaller and grows more slowly compared to that dominated by a mechanical thickening process. Thus, the water depth under the ice jam changes less. The depth of the scour hole around a pier develops similarly to that under an ice-covered flow condition, and the depth of the scour hole around a pier increases slowly. During the thickening process of the ice jam (dominated by a hydraulic thickening process) along the channel, the bottom of the ice jam appears as a wave-shaped accumulation. It is noticed that when the crest of such a wave migrates from upstream to the cross section of the pier, the flow depth under the ice jam decreases, and the flow velocity increases. Under such a circumstance (migration of a wave of ice accumulation), the scouring process at the pier becomes intense. Gradually, the development of the ice jam approaches an equilibrium condition, and

fewer ice particles entrained and accumulating under an ice cover/initial ice jam, the thickness of the initial ice jam is much smaller and grows more slowly compared to that dominated by a mechanical thickening process. Thus, the water depth under the ice jam changes less. The depth of the scour hole around a pier develops similarly to that under an ice-covered flow condition, and the depth of the scour hole around a pier increases slowly. During the thickening process of the ice jam (dominated by a hydraulic thickening process) along the channel, the bottom of the ice jam appears as a wave-shaped accumulation. It is noticed that when the crest of such a wave migrates from upstream to the cross section of the pier, the flow depth under the ice jam decreases, and the flow velocity increases. Under such a circumstance (migration of a wave of ice accumulation), the scouring process at the pier becomes intense. Gradually, the development of the ice jam approaches an equilibrium condition, and

**Figure 16.** Comparison of the development of a scour hole around a pier under open flow, ice cover and ice jam conditions in Case 2. **Figure 16.** Comparison of the development of a scour hole around a pier under open flow, ice cover and ice jam conditions in Case 2. **Figure 16.** Comparison of the development of a scour hole around a pier under open flow, ice cover and ice jam conditions in Case 2.

#### **4. Empirical Relationship between Scour Depth and Ice Jam Thickness 4. Empirical Relationship between Scour Depth and Ice Jam Thickness 4. Empirical Relationship between Scour Depth and Ice Jam Thickness**

Wang et al. (2021) performed experiments to study the local scour at a pier under an ice-jammed flow condition [3]. In their study, the relationship between the equilibrium thickness of ice jams (*Ti*) and the maximum scour depth of scour holes at the pier (*hs*) was derived as follows, Wang et al. (2021) performed experiments to study the local scour at a pier under an ice-jammed flow condition [3]. In their study, the relationship between the equilibrium thickness of ice jams (*Ti*) and the maximum scour depth of scour holes at the pier (*hs*) was derived as follows, Wang et al. (2021) performed experiments to study the local scour at a pier under an ice-jammed flow condition [3]. In their study, the relationship between the equilibrium thickness of ice jams (*T<sup>i</sup>* ) and the maximum scour depth of scour holes at the pier (*hs*) was derived as follows,

$$\frac{T\_i}{H} = 1.751 \left(\frac{h\_i}{H}\right) - 0.016 \quad \text{with } R^2 = 0.491 \tag{1}$$
 
$$\text{if the ratio of the maximum score depth of source holes to the flow}$$

 <sup>−</sup> 016.0751.1 *H* <sup>−</sup> 016.0751.1 <sup>=</sup> *<sup>H</sup> H* The dependence of the ratio of the maximum scour depth of scour holes to the flow depth (*H*) on the flow Froude number (*F<sup>r</sup>* = *V*/(*gH*) 0.5) can be expressed as the following:

$$\frac{\mu\_s}{H} = 2.014 \frac{V}{\sqrt{gH}} - 0.146 \quad \text{with } R^2 = 0.951 \tag{2}$$

where *V* = average flow velocity and *g* = gravitational acceleration.

In the above study [3], the impacts of both pier size and particle size of bed material on the maximum depth of scour holes at piers under an ice-jammed flow condition were not explored. The effect of the ice discharge rate on the maximum depth of scour holes at piers under an ice-jammed flow condition was not investigated. In addition, the shape of the piers was not considered.

Namaee and Sui [33–37] carried out experiments in a large outdoor flume to study local scour around two side-by-side piers under ice-covered flow conditions by using three nonuniform sands. Considering factors influencing the depth of scour holes around the side-by-side piers, such as flow Froude number, ice cover bottom roughness, bed roughness, pier diameter, and pier spacing, the following equation for describing the maximum scour depth of scour holes was obtained:

$$\frac{y\_{\text{max}}}{y\_0} = 5.96 \left(\frac{D\_{\text{50}}}{y\_0}\right)^{-0.070} \left(\frac{\text{G}}{\text{D}}\right)^{-0.256} \left(\frac{\text{n}\_{\text{i}}}{\text{n}\_{\text{b}}}\right)^{0.546} (\text{F}\_{\text{r}})^{1.67} \quad \text{with } \text{R}^2 = 0.900 \tag{3}$$

where *y*max = maximum scour depth, *D*<sup>50</sup> = median particle diameter, *y*<sup>0</sup> = approach flow depth, G = bridge spacing distance, *D* = pier diameter, *n*<sup>i</sup> = ice cover roughness, and *n*<sup>b</sup> = channel bed roughness.

In the above study [30], two types of model ice covers with different roughness coefficients were used. The evolution of an ice jam on the local scour process around piers was not studied.

Overall, the maximum depth of scour holes (*hs*) in the vicinity of the pier under an ice jam depends on the flow intensity, pier size (*D*), particle size of the sand bed (d*<sup>b</sup>* ), ice particle size (*d<sup>i</sup>* ), thickness of an ice jam (*T<sup>i</sup>* ), ice discharge rate (*Q<sup>i</sup>* ), roughness of the channel bed (*n<sup>b</sup>* ), and ice cover (*n<sup>i</sup>* ), respectively. The maximum depth of a scour hole (*hs*) can be expressed as the following relation:

$$h\_{\mathbb{S}} = f(T\_{\text{i}}, H, V, Q\_{\text{i}}, D, \underline{\text{g}}, d\_{\text{i}}, d\_{\text{b}}, n\_{\text{i}}, n\_{\text{b}}) \tag{4}$$

From the dimensional analysis, the relative maximum depth of a scour hole (*hs*/*H*) can be presented as follows,

$$\frac{h\_S}{H\_0} = f\left(\frac{T\_i}{H\_0}, \frac{V\_0}{\sqrt{gH}}, \frac{Q\_i}{Q\_w}, \frac{D}{H\_0}, \frac{d\_i}{d\_b}, \frac{n\_i}{n\_b}\right) \tag{5}$$

The dimensionless variable Froude number of the approaching flow is very important. In the current study, only one sand was used for bed material. For the model ice particle, only one size of model ice was used. Thus, only the dominant parameters affecting the maximum scouring depth around a pier under an ice jam were considered for developing an empirical relationship to describe the dependency of the maximum scour depth on the ice jam thickness, flow Froude number, and ice discharge rate. The ratio *ni*/*n<sup>b</sup>* was not considered in this study, assuming both *n<sup>i</sup>* and *n<sup>b</sup>* did not change much. In the present study, the pier coefficient was also considered since different pier shapes were used. Based on data collected from experiments, a regression analysis of the relative maximum scour depth (*hS*/*H*) against other dimensionless variables including the relative ice jam thickness (*Ti*/*H*), flow Froude number (*Fr*), grain size of bed material (*d*50/*D*), pier shape factor(*K<sup>ζ</sup>* ), and ice discharge rate (*Qi*/*Q*) was performed. The following equation for calculating the maximum depth of a scour hole at a pier under an ice jam was derived:

$$\frac{h\_6}{H\_0} = 10.928 \left(\frac{T\_i}{H\_0}\right)^{0.247} \left(\frac{V\_0}{\sqrt{gH\_0}}\right)^{1.483} \left(\frac{Q\_i}{Q\_w}\right)^{0.177} \left(\frac{d\_{\overline{g}0}}{D}\right)^{-0.416} \text{K}\_\zeta^{1.288} \quad \text{with } R^2 = 0.771 \tag{6}$$

where *g* = gravitational acceleration, *Q* = flow discharge, *H* = approaching flow depth, *V* = average approaching velocity, *d*<sup>50</sup> = Grain size of sand bed, *D* = pier size, and *K<sup>ζ</sup>* = pier shape factor.

Equation (3) implies that the relative maximum scour depth (*hS*/*H*) of scour holes around the pier primarily depends on the relative ice jam thickness (*Ti*/*H*), Froude number (*Fr* = *V*/(*gH*) 0.5) of the approaching flow, and ice discharge rate (*Qi*/*Q*). One can see from Equation (3) that the thicker the ice jam accumulation, the deeper the scour hole around the pier. The depth of a scour hole also increases with the flow Froude number. Interestingly, a larger ice discharge rate (*Qi*/*Q*) also results in a deeper scour hole. This is because a larger ice discharge rate (*Qi*/*Q*) leads to a thicker ice jam. On the other hand, the finer the bed material (*d*50/*D*), the larger and the deeper the scour hole around the pier. For a bridge pier whose cross section has a rectangular shape, the pier shape factor is more than that of a cylindrical shape. Thus, the maximum scour depth around the rectangular-shaped pier is more than that around a cylindrical pier. As shown in Figure 17, the results of the calculation using Equation (3) agree well with those of the laboratory measurements. The *hs* ' in Figure 17 is the calculated maximum depth of a scour hole. shaped pier is more than that around a cylindrical pier. As shown in Figure 11, the results of the calculation using Equation (3) agree well with those of the laboratory measurements. The *hs*' in Figure 17 is the calculated maximum depth of a scour hole.

288.1

ζ

416.0

*D*

− 

50

<sup>=</sup> with *R*2 = 0.771 (6)

*d*

177.0 483.1

*Q Q*

 

*w*

where *g* = gravitational acceleration, *Q* = flow discharge, *H* = approaching flow depth, *V* = average approaching velocity, *d*<sup>50</sup> = Grain size of sand bed, *D* = pier size, and *K<sup>ζ</sup>* = pier

Equation (3) implies that the relative maximum scour depth (*hS/H*) of scour holes around the pier primarily depends on the relative ice jam thickness (*Ti/H*), Froude number (*Fr* = *V*/(*gH*)*0.5*) of the approaching flow, and ice discharge rate (*Qi/Q*). One can see from Equation (3) that the thicker the ice jam accumulation, the deeper the scour hole around the pier. The depth of a scour hole also increases with the flow Froude number. Interestingly, a larger ice discharge rate (*Qi/Q*) also results in a deeper scour hole. This is because a larger ice discharge rate (*Qi/Q*) leads to a thicker ice jam. On the other hand, the finer the bed material (*d*50*/D*), the larger and the deeper the scour hole around the pier. For a bridge pier whose cross section has a rectangular shape, the pier shape factor is more than that of a cylindrical shape. Thus, the maximum scour depth around the rectangular-

0 0

*gH V*

928.10 *K*

 

*Water* **2022**, *14*, x FOR PEER REVIEW 17 of 20

247.0

*s i i*

   

0 0

*H h*

shape factor.

*H T*

 

### **5. Conclusions 5. Conclusions**

In the present study, experiments in a flume were carried out to study the influence of the evolution of ice jams on the local scour around bridge piers. The following results were obtained: In the present study, experiments in a flume were carried out to study the influence of the evolution of ice jams on the local scour around bridge piers. The following results were obtained:


**Author Contributions:** H.H.: Laboratory works, data curation, formal analysis, methodology, and writing—original draft preparation; J.W.: Conceptualization, laboratory supervision, methodology, funding acquisition, writing—original draft preparation; T.C.: Laboratory works, data curation, investigation; Z.H.: Laboratory works, data curation, investigation; J.S.: Conceptualization, formal analysis, methodology, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China (grant nos. 51879065). The authors are grateful for the financial support.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are available in the case that it is required.

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

### **References**


**Masoud Kazem <sup>1</sup> , Hossein Afzalimehr 1,\* and Jueyi Sui <sup>2</sup>**


**Abstract:** In presence of vegetation patches in a channel bed, different flow–morphology interactions in the river will result. The investigation of the nature and intensity of these structures is a crucial part of the research works of river engineering. In this experimental study, the characteristics of turbulence in the non-developed region downstream of a vegetation patch suffering from a gradual fade have been investigated. The changes in turbulent structure were tracked in sequential patterns by reducing the patch size. The model vegetation was selected carefully to simulate the aquatic vegetation patches in natural rivers. Velocity profile, TKE (Turbulent Kinetic Energy), turbulent power spectra and quadrant analysis have been used to investigate the behavior and intensity of the turbulent structures. The results of the velocity profile and TKE indicate that there are three different flow layers in the region downstream of the vegetation patch, including the wake layer, mixing layer and shear layer. When the vegetation patch is wide enough (*Dv*/*Dc* > 0.5, termed as the patch width ratio, where *D<sup>v</sup>* is the width of a vegetation patch and *D<sup>c</sup>* is the width of the channel), highly intermittent anisotropic turbulent events appear in the mixing layer at the depth of *z*/*H<sup>v</sup>* = 0.7~1.1 and distance of *x*/*Hv* = 8~12 (where *x* is streamwise distance from the patch edge, z is vertical distance from channel bed and *Hv* is the height of a vegetation patch). The results of quadrant analysis show that these structures are associated with the dominance of the outward interactions (Q1). Moreover, these structures accompany large coherent Reynolds shear stresses, anomalies in streamwise velocity, increases in the standard deviation of TKE and increases in intermittent Turbulent Kinetic Energy (TKE<sup>i</sup> ). The intensity and extents of these structures fade with the decrease in the size of a vegetation patch. On the other hand, as the size of the vegetation patch decreases, von Karman vortexes appear in the wake layer and form the dominant flow structures in the downstream region of a vegetation patch.

**Keywords:** coherent Reynolds shear stress; intermittent turbulence; mixing layer; quadrant analysis; vegetation patches

### **1. Introduction**

Submerged vegetation patches are fundamental participants of aquatic ecosystems, providing important habitats for fauna, including fishes [1–3]. As the reduction in the aquatic vegetation is a global concern, researchers have conducted studies to investigate the concept of a large-scale environmental issue [4–6]. However, while most of the current hydraulic-orientated research is focused on the flow characteristics in the canopy region of fully vegetated channels, finite dimensional vegetation patches with no fully developed region downstream of vegetation patches are very common in natural rivers. These sparse vegetation patches have a substantial effect on flow structures depending on the dimensions of the vegetation patch [7–9]. The effect of vegetation patch dimensions on flow structures has significant consequences for the ecosystem through various mechanisms. For example, an intermediate coverage (10–40%) of vegetation may promote high species richness

**Citation:** Kazem, M.; Afzalimehr, H.; Sui, J. Characteristics of Turbulence in the Downstream Region of a Vegetation Patch. *Water* **2021**, *13*, 3468. https://doi.org/10.3390/w13233468

Academic Editor: Achim A. Beylich

Received: 18 November 2021 Accepted: 1 December 2021 Published: 6 December 2021

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

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

due to the conditions provided by the optimal behavior of flow [10–12]. In addition, flow structures around a vegetation patch have a significant impact on the sediment transport, and consequently on the topography of the river bed [13–15]. Overall, the interaction between flow and vegetation in natural channels has attracted a lot of attention from hydraulic researchers. To mention only a few, the most recent research works have been conducted by researchers [16–18], in which various aspects of this interaction have been investigated.

Given the importance of the research topic, this present study aims to assess the effects of different sizes of vegetation patch on the flow characteristics in the downstream region of vegetation patches. The sizes of vegetation patches range from full coverage across the entire channel to a narrow vegetation patch located in the middle of the channel.

### **2. Flow Structures behind Vegetation Patches**

Compared to a rough bed without a vegetation patch, the presence of a vegetation patch in a channel results in changes in the flow structures and momentum transfer. Therefore, flow resistance and turbulence characteristics are expected to differ in the downstream regions of vegetation patches [19]. To assess the characteristics of turbulent flow, the following methodologies are commonly used: spectral analysis of velocity time series, determination of Reynolds stresses, and Turbulent Kinetic Energy (TKE), as well as quadrant analysis [20–22]. The analysis of the coherent Reynolds shear stress is another method that is used to determine coherent turbulent structures. However, this approach has hardly been used in hydraulic-oriented works due to the difficulties in calculating the phase velocity.

In order to make clear the terminologies used in the above-mentioned approaches, an introduction is presented in this section alongside the literature review. In the classic terminology of turbulences, for a sufficiently long duration of velocimetry (*T*), and for the velocity component of *u<sup>i</sup>* (*t*), the average velocity (*u<sup>i</sup>* ) is defined as:

$$\overline{u}\_{i} = \lim\_{T \to \infty} \frac{1}{T} \int\_{0}^{T} u\_{i}(t)dt \tag{1}$$

Considering *i = 1, 2, 3* as the *x*, *y* and *z* directions in the Cartesian coordinate system, *u*, *v* and *w* are defined as the averages of velocities in the *x, y*, and *z* directions, respectively. Consequently, the velocity fluctuations are defined by:

$$\begin{cases} u' = u(t) - \overline{u} \\ v' = v(t) - \overline{v} \\ w' = w(t) - \overline{w} \end{cases} \tag{2}$$

Then, the Turbulent Kinetic Energy (TKE) is also calculated as:

$$\text{TKE} = 0.5\overline{\left(u^{\prime 2} + \overline{v^{\prime 2}} + \overline{w^{\prime 2}}\right)}\tag{3}$$

In general, as was reported by researchers, a sharp increase in TKE is linked to the presence of vortexes [23–26]. In partially covered channels with emergent vegetation patches, where the flow is almost 2D, the flow structure in the downstream region of vegetation patches is dominated by the horizontal von Karman vortexes developed by instabilities in the trailing edges of the vegetation patch [27,28]. The distance between the trailing edge of the barrier and the formation of von Karman vortexes is known as *Lkv*. According to reported research based on TKE analysis, for a relatively dense vegetation patch with a diameter of *Dv*, in which the leakage velocity from the dense vegetation patch is negligible, *Lkv* ≈ 2.5*D<sup>v</sup>* [27,28]. On the other hand, the 3D flow structures in the downstream region of a submerged vegetation patch are more complicated. In addition to the constant presence of a vertical recirculation zone beyond the vegetation patch, the presence and magnitude of the von Karman vortexes in the horizontal plane depends on

the geometry of the vegetation patch [28–30]. For instance, the presence of von Karman vortexes at a depth of *Hv/H* > 0.55~0.7, where *H<sup>v</sup>* is the height of the vegetation patch and *H* is the flow depth [28,30]. On the other hand, results indicated that the dominant structure of a wide submerged path occurs in the vertical plane, as vortexes forming a vertical recirculation zone [28,29]. The distance between the trailing edge of the barrier and the forming of the vertical vortex is known as *Lkv*. For a solid submerged barrier, there is no gap between the trailing edge of the barrier and the forming vortex, and the length of wake behind the barrier is zero (*Lvv* = 0) [31]. However, for a porous vegetation patch, the velocity of flow passing through the vegetation patch may delay the formation of the recirculating region (*Lvv* 6= 0). The reported value of *Lvv* ranges between 1*H<sup>v</sup>* and 5*H<sup>v</sup>* depending on the porosity and diameter of the vegetation patch [28,29]. As reported by Liu et al. (2018), the value tends to be about 1*H<sup>v</sup>* in a relatively wide vegetation patch, and increases up to 5*H<sup>v</sup>* in vegetation patches with a low blockage ratio [28]. These results indicate that *Lvv* and *Lkv* are affected by the vegetation patch dimensions of *h* and *D*, respectively. In addition, the bleed velocity that penetrates into the wake region can also increase these values considerably [23,31].

In addition to the above-mentioned classic variables, coherent turbulent variables have become fundamental terms of flow analysis in recent years [32–34]. According to the triple decomposition approach, the instantaneous velocity can be written as [35]:

$$
\mu(\mathbf{x},t) = \overline{\mathbf{u}}(\mathbf{x}) + \tilde{u}\_{\mathbf{c}}(\mathbf{x},t) + u\_{r}(\mathbf{x},t) \tag{4}
$$

where *u*(*x*, *t*) is the instantaneous velocity, *u<sup>r</sup>* is the deflection of velocity (known as *u'* in the classic approach), and *u*(*x*) is the classic average of velocity during time period *T*. In addition to the classic averaging method, the periodic phase average can be defined as:

$$\langle \mu(\mathbf{x}, t) \rangle = \lim\_{N \to \infty} \frac{1}{T\_P} \sum\_{l}^{N} u\_l(\mathbf{x}, t + iT\_P) \tag{5}$$

where *T<sup>p</sup>* is the period of occurrence of a coherent structure and is equal to 1/*f<sup>d</sup> ,* where *f<sup>d</sup>* is the dominant frequency of the occurrence. Generally, when a coherent event occurs in a particular region of flow, *f<sup>d</sup>* will reach an obvious peak in the power density spectrum of the related time series in the same region. For a time series with a length of *T*, *N* is the number of cycles with a period of *Tp*, which can be calculated as *N = T/Tp*. Following these definitions, the coherent velocity deflection is:

$$
\widetilde{u}\_{\varepsilon}(\mathbf{x}, t) = \langle u(\mathbf{x}, t) \rangle - \overline{u}(\mathbf{x}) \tag{6}
$$

Considering the same approach for other components of velocity, the coherent Reynolds shear stress can be written as:

$$\pi\_{\mathfrak{c}} = -\rho \langle \tilde{u}\_{\mathfrak{c}} \tilde{w}\_{\mathfrak{c}} \rangle \tag{7}$$

The total coherent and non-coherent Reynolds shear stress is defined as follows:

$$
\pi\_r = -\rho \overline{\overline{\mu}\_r} \overline{\overline{w}\_r} \tag{8}
$$

Although it is very rare to implement *τ<sup>c</sup>* in studies of river hydraulics with the presence of vegetation patches, a few research works have been conducted to study the variation in *τ<sup>r</sup>* in river hydraulics in the presence of vegetation patches. For instance, in the downstream region of a submerged vegetation patch, there is a high gradient pattern of Reynolds shear stress in the vertical plane, which increases toward the surface of the flow [36]. This phenomenon is induced by the high-speed flow that passes above the submerged vegetation patch and causes a considerable velocity gradient when compared with the low velocity flow leaking through the vegetation patch. The region with a high velocity gradient is known as the shear layer. This region of flow has particular turbulence characteristics that will be studied in this paper.

As another method of analysis, the decomposition of bursting events (known as quadrant analysis in 2D space and octant analysis in 3D space) is also widely applied to determine the dominant turbulent events in the presence of both emerged and submerged vegetation patches [22,37–39]. The occurrence probability of event *k* is calculated as the normalized occurrence frequency, *f<sup>k</sup>* , for a particular class of events related to different classes of events [40]:

$$f\_k = \frac{n\_k}{N} \text{ with } N = \sum\_{1}^{4} n\_k \, k = 1, 2, 3, 4 \tag{9}$$

where *n<sup>k</sup>* is the number of events belonging to class *k*, and *N* is the total number of events.

For the quadrant analysis, a "hole" region has been proposed in the majority of previous studies. In a "hole" region, the event must be filtered and not be considered. In octant analysis, this threshold can be written as:

$$\left|u'(t)w'(t)\right| > \mathcal{C}\_H \left|\overline{u'w'}\right|\tag{10}$$

For the threshold parameter *CH*, the low-intensity events below a certain limit were omitted, which were scaled by the average of the velocity fluctuations. The high value of *C<sup>H</sup>* implies the selection of the strongest events, but the total number of instantaneous [*u'*(*t*) *w'*(*t*)] decreases so much that the contribution region in each quadrant becomes meaningless. For quadrant analysis, the threshold level *C<sup>H</sup>* = 1 is suggested for use to reach a good compromise between the clear identification of the events and the preservation of a number of instantaneous events of each class [40,41]. In this paper, *C<sup>H</sup>* is set as 1. The classes of events are also defined as:


To date, most of the reported research has applied quadrant analysis to study coherent flow structures within or above a vegetation patch. However, there are a few published studies in which quadrant analysis has been conducted to assess the coherent flow structures in the downstream region of a vegetation patch. Okamoto and Nezu (2013) point out that in the region immediately behind the submerged vegetation patch, there is a small margin of ejection dominancy around the top edge [7]. Similarly, based on data collected from experiments performed on dryland vegetation in a wind tunnel, Mayaud et al. (2016) revealed that there were elevated frequencies of Q2 (ejection) and Q4 (sweep) events in the immediate toe of the vegetation patch [42]. In contrast, it is reported the dominance of outward and inward interactions in the shear layer induced by the flow passing above the vegetation patch, which is a notable characteristic of the flow in the downstream region of a submerged vegetation patch in a more distant region [38]. It must be mentioned that the relation between quadrant occurrences and velocity structures is an interesting research topic that has received attention from researchers recently. Wang et al. (2019) established interconnections between the classic definition of vortex groups and quadrant occurrences, which are used in the next section of this paper [43].

In conclusion, the introduced techniques have been successfully used in the hydrodynamics analysis in channels with the presence of vegetation patches. However, the reported studies on the hydrodynamics of channels with submerged vegetation patches are limited to those with the flow either above the patches or immediately downstream of the patches. For example, the flow structures was traced up to a distance of 8*H<sup>v</sup>* from the toe of the patch towards the downstream [36]. Liu et al. (2018) studied the flow structures to a distance of just 5~6*D<sup>v</sup>* from the patch toe in the downstream direction [28]. The quadrant analysis has been conducted in the field; here, the distance from the patch toe was even more limited, to ~3*D<sup>v</sup>* [38]. Although this distance from the patch toe is enough to address the characteristics of patch-induced turbulence in the near field, it is

not enough to cover the more distant events of turbulence. For example, the investigation of the characteristics of the shear layer generated above the patch requires an extended range of velocimetry in the downstream direction. In addition, there is a need to combine the above-mentioned methods to evaluate the turbulent characteristics of flow behind the submerged vegetation patches. This requires a particular emphasis on the implementation of *τ<sup>c</sup>* in the turbulence analysis, which has hardly been studied in recent works. In order to fill the above-mentioned gaps in the previous research, in the present study, a combined approach to turbulent analysis is used to investigate the turbulent characteristics of flow in the region downstream of a vegetation patch with an extended range up to 17*Hv*.

### **3. Experimental Setup**

The experiments were conducted using a glass flume at the Laboratory of Hydraulics of the Iran University of Science and Technology. The flume is 14 m long, 0.9 m wide and 0.6 m deep. The discharge was controlled by an electromagnetic flow meter installed at the entrance of the flume and was set for 31 L/s in the present study. The water level in the flume was adjusted by a tailgate located at the end of the flume and was set for a depth of 18.5 ± 0.3 cm in the present study. The distance between the flume entrance and vegetated flume zone was 6 m to ensure a fully developed flow in the region upstream of the vegetation patch.

The velocity measurements were conducted after the flow reached the steady state condition. The velocity profiles were measured using an acoustic Doppler velocimeter (ADV), placed at the centerline of each row of vegetation patch. There were 19–26 measuring points along each vertical line for velocity measurements, and the vertical distance between two adjacent measuring points was 4–10 mm. The sampling frequency and the measuring time of the ADV were 200 Hz and 120 s, respectively, resulting in 24,000 instantaneous velocities for point measurement.

To date, most previous studies have employed artificial simple rods of regular shapes to simulate natural vegetation patches. As real vegetation is flexible and irregular, this method might not represent the nature of vegetation's behavior [36]. On the other hand, the experiments showed that natural vegetation would lose its stiffness and develop a long-lasting curvature in the flow direction after couple of days. Therefore, a well-shaped synthetic plant was used to represent the natural vegetation, which was selected on the basis of a real world sample of patches in a gravel bed river.

Each model plant had three branches. Each branch had 12 leaves, and the diameter of the branch trunk was approximately 3 mm, as shown in Figure 1. The average height of the vegetation patch was 105 ± 5 mm, and the lateral and longitudinal spread widths of the leaves were approximately 9–19 and 11–22 mm, respectively. The model plants had a certain degree of flexibility and could swing in a flowing current similar to vegetation in a natural river. The vegetation patch was attached to a perforated board in a staggered arrangement. Four different layouts of vegetation patches were used to simulate both fully covered and non-fully covered channel beds. The dimensions of these vegetation patches, namely, length (*Lv*) × width (*Dv*), in this study were 120 × 90 cm, 90 × 60 cm, 60 × 45 cm and 40 × 30 cm, respectively. Table 1 summarizes the experimental runs. The material on the flume bed was a mixture of natural gravel similar to that in a natural gravel bed river (the Marbor River, Zagros Mountains region, Iran). The equivalent particle diameter, which 90% of the total particles were smaller than (d90), was 18.8 mm, and d<sup>50</sup> = 14 mm. The layout of the experimental device and employed materials are shown in Figure 1.

**Table 1.** Experimental parameters for the cases discussed in this paper.

**n/m2 (Number of Veg. per Square Meter)**

31 L/s 611.1 120 90 10 31 L/s 611.1 90 60 10 31 L/s 611.1 60 45 10 31 L/s 611.1 40 30 10 **No. Veg.** 31 L/s - - - -

**Lv (Length of Patch, cm)** 

**Dv (Width of Patch, cm)** 

**Hv (Height of Patch, cm)** 

**Case Q (Discharge,** 

**L/s)** 

**Figure 1.** (**A**) The layout of the experimental device, (**B**) a sample of vegetation patches (120 × 90 cm, *DV/DC* = 1), (**C**) sample of vegetation patch in a natural gravel bed river (Marbor River, Iran), (**D**) a sample of bed materials, (**E**) a single synthetic plant used to simulate patch; (**F**) a sketch of the experimental setup (not to scale). **Figure 1.** (**A**) The layout of the experimental device, (**B**) a sample of vegetation patches (120 × 90 cm, *Dv/D<sup>c</sup>* = 1), (**C**) sample of vegetation patch in a natural gravel bed river (Marbor River, Iran), (**D**) a sample of bed materials, (**E**) a single synthetic plant used to simulate patch; (**F**) a sketch of the experimental setup (not to scale).



### structures in the downstream region of a vegetation patch is provided. **4. Results and Discussion**

The results are presented as follows: First, both the velocity profile and the TKE have been investigated to determine the flow layers that formed behind the vegetation patch. Then, spectral analysis and coherent Reynolds shear stress analysis have been carried out to determine the nature of turbulence and coherent Reynolds shear stresses in different flow layers. Consequently, the results of the quadrant analysis of bursting events are provided to compare and evaluate the coherent occurrences in the framework of the time domain. Then, by considering the particular characteristics of turbulence in the mixing layer, the temporal characteristics of turbulence in the mixing layer are investigated. Finally, by

nally, by considering the results of previous sections, the transformation of coherent

considering the results of previous sections, the transformation of coherent structures in the downstream region of a vegetation patch is provided.

### *4.1. Velocity Profile and TKE*

For all experimental runs, three layers of flow can be observed based on the inflection points of the velocity profiles (Figure 2). The associated regions of these layers are known as the wake zone, the mixing layer and the Log–Law shear zone. These three regions resemble the results of trough patch velocimetry reported by other researchers [44–47]. For instance, in case 1 and for a distance from the vegetation edge *x/H<sup>v</sup>* = 8, the wake zone formed at the water depth of *z/H<sup>v</sup>* = 0~0.7 (*z/H* = 0~0.56), the mixing layer formed at the water depth of *z/H<sup>v</sup>* = 0.7~1.1 (*z/H* = 0.56~0.62) and the Log–Law shear zone formed at the water depth of *z/H<sup>v</sup>* > 1.1 (*z/H* > 0.62). The extent and thickness of the mixing layer reduced as the vegetation patch decayed. In this layer, a severe deflection is detectable in the velocity profile, which resembles an adverse pressure gradient. For case 1, this anomalous phenomenon covered a distance (from the vegetation edge) between 8 < *x/H<sup>v</sup>* and *x/H<sup>v</sup>* ≥ 17. In case 2, this phenomenon covered a distance of 3 < *x/H<sup>v</sup>* < 8, and it was only visible around *x/H<sup>v</sup>* = 12 in case 3. However, this phenomenon was absent in case 4. The wake zone was formed in the lower region beyond the patch. For case 1, the wake zone occurred at a distance from the vegetation edge of *x/H<sup>v</sup>* = 1, with the upper limit of this layer at *z/H* = 0.45 or *z/H<sup>v</sup>* = 0.8, which increased up to *z/H<sup>v</sup>* = 1 for the smallest patch (case 4). In this layer, the velocity profiles tended to adopt a Log–Law shape beyond the distance of *x/H<sup>v</sup>* = 12 for all cases. For a channel bed partially covered by a vegetation patch, the thickness of the wake layer varied along the flow direction, and increased slightly in such a way that it reached the lower boundary of the Log–Law shear region. Consequently, the mixing layer became narrower and even disappeared in the presence of smaller patches. However, for the case of a channel bed fully covered with a vegetation patch (case 1), the thickness of the mixing layer increased considerably at a distance (from the vegetation edge) of *x/H<sup>v</sup>* = 8, and seemed to be effectively present at a distance of *x/H<sup>v</sup>* >> 17.

While the vertical vortexes that normally generate behind a barrier have attracted a lot of attention from researchers, the development of horizontal vortexes (the von Karman Vortex Street) needs to be investigated too. Generally, the formation of von Karman vortexes is associated with the occurrence of peaks in the TKE values [23]. Thus, for each case, the TKE values have been calculated for the relative flow depths of *z/H<sup>v</sup>* = 0.5 and *z/H<sup>v</sup>* = 1, and also for the near-bed region behind the vegetation patch along the centerline of the flume (Figure 3). The variation in the streamwise velocity along the flow direction is also shown in this figure. In case 1, there was only one peak in the graph of TKE/(*U*0) <sup>2</sup> around *x/D<sup>v</sup>* = 0.9 and *x/H<sup>v</sup>* = 7.5~8, which is associated with severe deflection in the velocity profile. However, no flow escaped from the sides of the vegetation patch (since the channel bed was fully covered with vegetation). This peak in TKE/(*U*0) 2 clearly indicates the presence of a vertical vortex behind the vegetation patch and near the water's surface. While a sharp increase in TKE occured after a distance of *x/H<sup>v</sup>* = 4.5~5, the length of the wake behind the patch can be considered as *Lvv* = 5*Hv*. Importantly, the anomalies in the mixing layer appeared in the same point as was reported [28] for a length of wake of *Lvv* = (3.5~5)*H<sup>v</sup>* with a wide submerged vegetation patch. However, Folkard (2005) reported a slightly lower range for a wide but highly submerged vegetation patch [29]. Overall, these results are comparable to those of Liu et al. (2018) for a wide and partially vegetation-covered bed with slight submergence [28].

*4.1. Velocity Profile and TKE* 

>> 17.

For all experimental runs, three layers of flow can be observed based on the inflection points of the velocity profiles (Figure 2). The associated regions of these layers are known as the wake zone, the mixing layer and the Log–Law shear zone. These three regions resemble the results of trough patch velocimetry reported by other researchers [44– 47]. For instance, in case 1 and for a distance from the vegetation edge *x/Hv* = 8, the wake zone formed at the water depth of *z/Hv* = 0~0.7 (*z/H* = 0~0.56), the mixing layer formed at the water depth of *z/Hv* = 0.7~1.1 (*z/H* = 0.56~0.62) and the Log–Law shear zone formed at the water depth of *z/Hv* > 1.1 (*z/H* > 0.62). The extent and thickness of the mixing layer reduced as the vegetation patch decayed. In this layer, a severe deflection is detectable in the velocity profile, which resembles an adverse pressure gradient. For case 1, this anomalous phenomenon covered a distance (from the vegetation edge) between 8 < *x/Hv* and *x/Hv* ≥ 17. In case 2, this phenomenon covered a distance of 3 < *x/Hv* < 8, and it was only visible around *x/Hv* = 12 in case 3. However, this phenomenon was absent in case 4. The wake zone was formed in the lower region beyond the patch. For case 1, the wake zone occurred at a distance from the vegetation edge of *x/Hv* = 1, with the upper limit of this layer at *z/H* = 0.45 or *z/Hv* = 0.8, which increased up to *z/Hv* = 1 for the smallest patch (case 4). In this layer, the velocity profiles tended to adopt a Log–Law shape beyond the distance of *x/Hv* = 12 for all cases. For a channel bed partially covered by a vegetation patch, the thickness of the wake layer varied along the flow direction, and increased slightly in such a way that it reached the lower boundary of the Log–Law shear region. Consequently, the mixing layer became narrower and even disappeared in the presence of smaller patches. However, for the case of a channel bed fully covered with a vegetation patch (case 1), the thickness of the mixing layer increased considerably at a distance (from the vegetation edge) of *x/Hv* = 8, and seemed to be effectively present at a distance of *x/Hv*

**Figure 2.** Streamwise velocity profiles along the centerline of the flume; *x/Hv* is the normalized distance from the downstream edge of the vegetation patch, and *H* is the depth of flow. The **Figure 2.** Streamwise velocity profiles along the centerline of the flume; *x/H<sup>v</sup>* is the normalized distance from the downstream edge of the vegetation patch, and *H* is the depth of flow. The boundaries between zones affected by the wake and over the canopy flow are shown as red lines. The velocity of the Log–Law zone is shown as dashed and is not to scale.

For experiment cases 2, 3 and 4, the location of the minimum velocity was shifted towards the flow direction. However, since the patch width ratio of *Dv/D<sup>c</sup>* was relatively high for case 2, *Lvv* was only shifted to the near-bed region. For cases 3 and 4, the shift of wake was also detectable at depths of *z/H<sup>v</sup>* = 0.5 and *z/H<sup>v</sup>* = 1. In these cases, the shift extent of the upper layers exceeded that of the near-bed zone. For all four cases, the minimum velocity was negative in the near-bed region, implying the occurrence of a vertical vortex. Considering the velocity values in the upper zone, it was observed that the center of the vortex (or rotation) moved upward if *Dv/D<sup>c</sup>* decreased. In case 1, the center of vortex was located in the near-bed region and *z/H<sup>v</sup>* = 0.5, while it was located between *z/H<sup>v</sup>* = 0.5 and *z/H<sup>v</sup>* = 1 in case 4 (smallest vegetation patch). This elevated center of rotation in case 4 pushed the larger part of flow down into the lower zones. Subsequently, the velocity pattern in the middle region became similar to that of the near-bed region (one can observe this by comparing the minimum velocity located between the near-bed region and *z/H<sup>v</sup>* = 0.5). In addition, the larger radius of the rotation in the channel with a smaller patch width ratio of *Dv/D<sup>c</sup>* caused the minimum streamwise velocity to be in the near-bed region rather than in the mixing layer. The reason for the upward movement of the rotation center in the case of a smaller patch in the channel may be the higher through-patch velocity. Regarding the larger values of *Lvv* due to the through-patch flow, this statement is consistent with the findings of other researchers [27,28,31]. In addition, for experiment cases 2, 3 and 4, the flow from the sides of the vegetation patches created the second peak in the TKE curve at the depth of *z/H<sup>v</sup>* = 1. However, for cases 2 and 3, while the *Dv/D<sup>c</sup>* values were still relatively high, the second peak in the TKE curve did not occur in the middle and near-bed region. Additionally, the location of the single peak in the TKE curve in the lower zones was shifted forward.

**Figure 3.** Normalized TKE and streamwise velocity behind 4 different patch layouts near channel bed (z/Hv = 0), at the middle height of patch (z/Hv = 0.5) and top of patch (z/Hv = 1). **Figure 3.** Normalized TKE and streamwise velocity behind 4 different patch layouts near channel bed (*z/Hv* = 0), at the middle height of patch (*z/Hv* = 0.5) and top of patch (*z/Hv* = 1).

For experiment cases 2, 3 and 4, the location of the minimum velocity was shifted towards the flow direction. However, since the patch width ratio of *Dv/Dc* was relatively high for case 2, *LVV* was only shifted to the near-bed region. For cases 3 and 4, the shift of wake was also detectable at depths of *z/Hv* = 0.5 and *z/Hv* = 1. In these cases, the shift extent of the upper layers exceeded that of the near-bed zone. For all four cases, the minimum velocity was negative in the near-bed region, implying the occurrence of a vertical Overall, the results of both the TKE and velocity analyses show that the three layers (wake, mixing and shear) were present in all cases of vegetation patches. However, the wake and shear layers were considerably affected by the flow passing from the sides of the vegetation patches, and as a consequence, different flow structures and associated length scales emerged, dependent on the upward movement of the center of the vertical rotation beyond the patch.

vortex. Considering the velocity values in the upper zone, it was observed that the center

### *4.2. Spectral Analysis and Coherent Reynolds Shear Stress*

The dominant frequency of the streamwise Reynolds shear stress (*fd<sup>u</sup>* 0*w*0 ) was the peak of the power spectrum diagram calculated via *u* 0*w* 0 fluctuations (*Suw*). Therefore, *Suw* is calculated at different points for all experimental cases (Figure 4), and the results are used to calculate the coherent Reynolds shear stresses (−h*u*e*cw*e*c*i). The results of the spectral analysis can be placed in the following three categories:


tropic turbulence.

zones and the center of rotation developed in the upper zone. Such a weak vortex cannot produce the strong vertical vortexes required to alter the isotropic turbulence. A Matlab@ code was used to calculate the phase velocity and its deflection for the dominant frequency at each point. The coherent vertical momentum transfer triggered by coherent vertical vortexes continued to appear in both the near-canopy and near-bed regions in cases 1 to 3. However, near-bed coherent shear stress did not occur in case 4. It can be inferred that the flow regime in the near bed region of case 4 was dominated by the flow through the vegetation patch and the wall effect of the rough bed. The presence of the vegetation patch had no coherent effect on this region. In addition, for cases 1, 2 and 3, the majority of coherent occurrences were accompanied by a high spatial gradient of the Reynolds shear stress, which is associated with the boundary layer separation zone and large scale vortexes—a characteristic that was described by Lian (1990) [48]. However, in case 1, the occurrence of coherent Reynolds shear stresses was associated with very low Reynolds shear stresses with a low spatial gradient. Consequently, it can be referred that the coherent Reynolds shear stresses can be classified into two categories according to the origins of the coherence. The first and most prevalent category is the strong coherent Reynolds stress, which was prevalent in the mixing layers of cases 1, 2 and 3. The second category is the weak coherent Reynolds shear stress that occurred in the wake layer. Based on the results of the spectra analysis, the strong coherent shear stresses were associated with peaks of the anisotropic turbulent spectra. In contrast, the weak form of the coherent shear stress occurred in the presence of the peak iso-

**Figure 4.** Samples of the power spectra behind the patch showing different classes of turbulence: (A) isotropic–homogenious turbulent flow without dominant frequency on *Suw*, (**B**) isotropic– homogenious turbulent flow with dominant frequency on Suw, observed below the height of the patch, (**C**) anisotropic turbulent flow with dominant frequency on Suw, observed alongside the height of patch, (**D**) Suw of the same anisotropic sample. **Figure 4.** Samples of the power spectra behind the patch showing different classes of turbulence: (**A**) isotropic–homogenious turbulent flow without dominant frequency on *Suw*, (**B**) isotropic– homogenious turbulent flow with dominant frequency on *Suw*, observed below the height of the patch, (**C**) anisotropic turbulent flow with dominant frequency on *Suw*, observed alongside the height of patch, (**D**) *Suw* of the same anisotropic sample. *Water* **2021**, *13*, x FOR PEER REVIEW 12 of 23

**Figure 5.** Values of −′′ തതതതതത along the centerline of the channel beyond the vegetation patch (colored figures) beside values of / (graysclae). Note that the blank area in / is where no dominant frequency was detectable in Suw. The end of the patch is located at x = 800. **Figure 5.** Values of −*u*0*w*0 along the centerline of the channel beyond the vegetation patch (colored figures) beside values of *τc*/*τ<sup>r</sup>* (graysclae). Note that the blank area in *τc*/*τ<sup>r</sup>* is where no dominant frequency was detectable in *Suw*. The end of the patch is located at *x* = 800.

On the other hand, for the patch width ratio of Dv/Dc = 0.5, a new type of coherent structure formed in the wake zone, and extended as the patch size reduced. With the patch width ratio of Dv/Dc = 0.33 (case 4), the greatest part of wake zone was covered by On the other hand, for the patch width ratio of *Dv/D<sup>c</sup>* = 0.5, a new type of coherent structure formed in the wake zone, and extended as the patch size reduced. With the patch width ratio of *Dv/D<sup>c</sup>* = 0.33 (case 4), the greatest part of wake zone was covered

this structure. While, in this region, the turbulent spectra of instantaneous velocity re-

from the patch sides could form strong von Karman vortexes in the horizontal plane. As reported by Siddique et al. (2008) and Liu et al. (2018), von Karman vortexes appeared in partially covered patches with *Hv/H* > 0.55~0.7 [28,30]. Consequently, for these regions, the frequency of vortex shedding was calculated as 0.07~0.11Hz, which is equivalent to an average Strouhal number of about 0.25. These results are comparable with those re-

The occurrence probabilities of the quadrant classes are shown in Figure 6. One can see from Figure 6 that in the presence of smaller patches, the ejection-dominated zone was shifted to a higher elevation at a shorter distance from the patch. This phenomenon was accompanied with lower Reynolds shear stress in the near-bed region. It was also associated with a marginally thicker wake layer in the presence of smaller patches, which

ported by other researchers [27,28].

*4.3. Quadrant Analysis of Bursting Events* 

was described in the previous section.

by this structure. While, in this region, the turbulent spectra of instantaneous velocity resembled isotropic turbulence, the effect of the interaction between the mixing layer and the shear layer seems to be insignificant. On the other hand, the presence of strong flow from the patch sides could form strong von Karman vortexes in the horizontal plane. As reported by Siddique et al. (2008) and Liu et al. (2018), von Karman vortexes appeared in partially covered patches with *Hv/H* > 0.55~0.7 [28,30]. Consequently, for these regions, the frequency of vortex shedding was calculated as 0.07~0.11 Hz, which is equivalent to an average Strouhal number of about 0.25. These results are comparable with those reported by other researchers [27,28].

### *4.3. Quadrant Analysis of Bursting Events*

The occurrence probabilities of the quadrant classes are shown in Figure 6. One can see from Figure 6 that in the presence of smaller patches, the ejection-dominated zone was shifted to a higher elevation at a shorter distance from the patch. This phenomenon was accompanied with lower Reynolds shear stress in the near-bed region. It was also associated with a marginally thicker wake layer in the presence of smaller patches, which was described in the previous section. *Water* **2021**, *13*, x FOR PEER REVIEW 13 of 23

**Figure 6.** Probability of occurrence of quadrant classes behind the patch; the end of the patch is located at x = 800. The height of the patches is 105 ± 5 mm, and the values are shown up to level of the maximum domain of ADV capability (~140 mm). **Figure 6.** Probability of occurrence of quadrant classes behind the patch; the end of the patch is located at *x* = 800. The height of the patches is 105 ± 5 mm, and the values are shown up to level of the maximum domain of ADV capability (~140 mm).

In contrast, the sweep-dominated zone, observed above the top of the patch in case 1, faded with the decrease in the patch dimensions. As regards the Reynolds shear stresses of this region, the occurrences of the sweep were associated with the high Reyn-

As strong coherent shear stresses were prevalent in the boundary between the mixing layer and the shear Log–Law layer, the quadrant analysis of the bursting events of this zone are outlined in Figure 7. The propinquity between the outward-dominated regions and regions with coherent Reynolds shear stresses is another notable result of the quadrant analysis. For all vegetation patches, the coherent shear stress in the mixing layer with a frequency of fୢ୳ᇱ୵ᇱ appeared in exactly the same region in which the outward class of quadrant occurrences were dominant, where *fk=1* > 0.25. This was a one-way relationship, and there were some points with *fk=1* > 0.25 at which no coherent shear stress was observable. To assess this important outcome, for point and quadrant class , two

A(x, i) = ൜1 if f୩ୀ୧ = maxሼf୩ୀଵ..ସሽ

0 if f୩ୀ୧ ≠ maxሼf୩ୀଵ..ସሽ

(11)

pattern was reported in the downstream region of dryland vegetation patches facing

conditional functions can be defined, as below:

wind flow [42].

In contrast, the sweep-dominated zone, observed above the top of the patch in case 1, faded with the decrease in the patch dimensions. As regards the Reynolds shear stresses of this region, the occurrences of the sweep were associated with the high Reynolds shear stress of the Log–Law layer in the downstream region of the patches. The same pattern was reported in the downstream region of dryland vegetation patches facing wind flow [42]. *Water* **2021**, *13*, x FOR PEER REVIEW 14 of 23 B(x) = ቊ1 if fୢ(౫ᇲ౭ᇲ) ≠ 0 0 if f୩ୀ୧ = 0 (12)

As strong coherent shear stresses were prevalent in the boundary between the mixing layer and the shear Log–Law layer, the quadrant analysis of the bursting events of this zone are outlined in Figure 7. The propinquity between the outward-dominated regions and regions with coherent Reynolds shear stresses is another notable result of the quadrant analysis. For all vegetation patches, the coherent shear stress in the mixing layer with a frequency of *fd<sup>u</sup>* 0*w*0 appeared in exactly the same region in which the outward class of quadrant occurrences were dominant, where *f <sup>k</sup>*=1 > 0.25. This was a one-way relationship, and there were some points with *f <sup>k</sup>*=1 > 0.25 at which no coherent shear stress was observable. To assess this important outcome, for point *x* and quadrant class *i*, two conditional functions can be defined, as below: Subsequently, ρ, (the Pearson correlation coefficient between A(x, i) and (x)) can be used to determine which class of quadrant occurrences is well-matched with the coherent Reynolds shear stress. For a sample including 278 points for all four cases of patch setup, ρ, is 0.68 for Q1 (outward interaction). In comparison, this value is −0.45, 0.26 and −0.38 for Q2, Q3 and Q4, respectively. This result confirms that the outward interactions are the representative quadrant class of the coherent Reynolds shear stresses in the downstream region of vegetation patches, particularly for strong coherent occurrences. Figure 8 illustrates the values of f୩ and fୢ୳ᇱ୵ᇱ for the 278 points beyond the vegetation patches that were investigated in this research. According to the definition of the out-

$$\mathbf{A}(\mathbf{x},i) = \begin{cases} 1 \text{ if } f\_{k=i} = \max\{f\_{k=1..A}\} \\ 0 \text{ if } f\_{k=i} \neq \max\{f\_{k=1..A}\} \end{cases} \tag{11}$$

$$\mathsf{B}(\mathsf{x}) = \left\{ \begin{array}{l} 1 \text{ if } f\_{d\_{(u'w')}} \neq 0 \\ 0 \text{ if } f\_{k=i} = 0 \end{array} \right. \tag{12}$$

**Figure 7.** Conceptual illustration of the quadratic bursting events in the boundary between the mixing layer and the shear layer. **Figure 7.** Conceptual illustration of the quadratic bursting events in the boundary between the mixing layer and the shear layer.

Subsequently, ρA,B (the Pearson correlation coefficient between A(*x*, *i*) and (*x*)) can be used to determine which class of quadrant occurrences is well-matched with the coherent Reynolds shear stress. For a sample including 278 points for all four cases of patch setup, ρA,B is 0.68 for Q1 (outward interaction). In comparison, this value is −0.45, 0.26 and −0.38 for Q2, Q3 and Q4, respectively. This result confirms that the outward interactions are the representative quadrant class of the coherent Reynolds shear stresses in the downstream region of vegetation patches, particularly for strong coherent occurrences. Figure 8 illustrates the values of *f<sup>k</sup>* and *fdu* <sup>0</sup>*w*<sup>0</sup> for the 278 points beyond the vegetation patches that were investigated in this research. According to the definition of the outward quadrant, the upward momentum flux within the mixing layer toward the shear layer must be considered

(see Figure 7). This upward flux is also detectable in the velocity profiles (Figure 2), where anomalies in the mixing layer are exacerbated; the streamwise velocity was reduced in the mixing layer and increased in the bottom part of the shear layer. This phenomenon illustrates a particular behavior in the velocity time series, which will be discussed in the next section. *Water* **2021**, *13*, x FOR PEER REVIEW 15 of 23

**Figure 8.** Values of (1) and fୢ୳ᇱ୵ᇱ for 278 points in the downstream region of the vegetation patches. **Figure 8.** Values of *f<sup>k</sup>* (Q1) and *fd<sup>u</sup>* 0*w*0 for 278 points in the downstream region of the vegetation patches.

### *4.4. Temporal Characteristics of Turbulence in the Mixing Layer 4.4. Temporal Characteristics of Turbulence in the Mixing Layer*

As was described in the previous sections, the mixing layer beyond the vegetation patch is prone to anisotropic turbulence, characterized by outward interaction events and the existence of coherent Reynolds stresses. The deviation from an ideal isotropic turbulence originated from the turbulent intermittency. The turbulent intermittency was characterized by a preference for turbulence with large velocity gradients, which is reflected in the strongly non-Gaussian tails of the probability density functions of velocity differences. These tails are determined by extreme events [49]. Intermittency is also defined as an abnormality in turbulent flow initiated by the interaction between turbulent regions or the interaction between a turbulent region and a vicinal laminar behavior region. The second mechanism is the common source of intermittency in shear flows. Moreover, the occurrence of a highly intermittent region in the near-boundary regions was associated with the prevalence of coherent vortices [50]. The mechanism generates quasi-laminar intervals in the velocity time series of a turbulent flow, which postpones the non-viscous intermediate sub-range in the energy cascade. Thus, there was a slight slope in the energy spectra. Consequently, these intervals persist up to arbitrarily large Reynolds numbers through the energy cascade sequences, which produce extreme events As was described in the previous sections, the mixing layer beyond the vegetation patch is prone to anisotropic turbulence, characterized by outward interaction events and the existence of coherent Reynolds stresses. The deviation from an ideal isotropic turbulence originated from the turbulent intermittency. The turbulent intermittency was characterized by a preference for turbulence with large velocity gradients, which is reflected in the strongly non-Gaussian tails of the probability density functions of velocity differences. These tails are determined by extreme events [49]. Intermittency is also defined as an abnormality in turbulent flow initiated by the interaction between turbulent regions or the interaction between a turbulent region and a vicinal laminar behavior region. The second mechanism is the common source of intermittency in shear flows. Moreover, the occurrence of a highly intermittent region in the near-boundary regions was associated with the prevalence of coherent vortices [50]. The mechanism generates quasi-laminar intervals in the velocity time series of a turbulent flow, which postpones the non-viscous intermediate sub-range in the energy cascade. Thus, there was a slight slope in the energy spectra. Consequently, these intervals persist up to arbitrarily large Reynolds numbers through the energy cascade sequences, which produce extreme events in velocity time series [51]. The determination of a particular intermittency index is a challenging issue and requires the spatial analysis of the velocity time series measured at different locations at the same time. Clearly, there is a need for appropriate sensors, such as PIV and hot wires.

in velocity time series [51]. The determination of a particular intermittency index is a challenging issue and requires the spatial analysis of the velocity time series measured at different locations at the same time. Clearly, there is a need for appropriate sensors, such as PIV and hot wires. However, temporal analyses of velocity time series can also provide a proper description of the quantity and intensity of intermittent occurrences. As the intermittency is linked with the non-Gaussian tails of the probability density functions of velocity differences, the velocity time series for the mixing layer were subjected to a normality test. The traditional

However, temporal analyses of velocity time series can also provide a proper description of the quantity and intensity of intermittent occurrences. As the intermittency is

ferences, the velocity time series for the mixing layer were subjected to a normality test. The traditional empirical rule of normality (or a simple normality test) states that, in a normally distributed data set, 68%, 95%, and 99.7% of the values lie within one, two, and three standard deviations (σ) of the mean (μ), respectively. Most of the points in the wake area could barely pass these restrictions; however, the normality test of Kolmogorov– Smirnov identified almost all points as non-normally distributed data sets. Regarding the methodology and application of the Kolmogorov–Smirnov normality test, please refer to

empirical rule of normality (or a simple normality test) states that, in a normally distributed data set, 68%, 95%, and 99.7% of the values lie within one, two, and three standard deviations (σ) of the mean (µ), respectively. Most of the points in the wake area could barely pass these restrictions; however, the normality test of Kolmogorov–Smirnov identified almost all points as non-normally distributed data sets. Regarding the methodology and application of the Kolmogorov–Smirnov normality test, please refer to Massey (1951) [52]. The results of the simple and Kolmogorov–Smirnov normality tests for case 1 are shown in Table 2. As the Kolmogorov–Smirnov normality test is very sensitive to non-Gaussian deviations, it even identifies semi-normal velocity time series for the wake layer as non-normal distributions. However, according to the concept of the Kolmogorov–Smirnov normality test, the intensity of non-normality can be observed on the basis of the deviation from the normal cumulative distribution. In this way, an error index can be defined as:

$$Error\_{\rm CDF} = 1 - \frac{CDF\_{\rm E}}{CDF\_{\rm N}} \tag{13}$$

where *CDF<sup>E</sup>* refers to the cumulative distribution of the empirical velocity time series, and *CDF<sup>N</sup>* is the cumulative distribution of a normal data set with the same µ and σ as the empirical data set. Figure 9 shows the variation in *ErrorCDF* for different depths at a distance of *x/H<sup>v</sup>* = 8, where the mixing layer is well developed. One can see that the values of *ErrorCDF* are considerably higher at a depth of 0.8 < *z/H<sup>v</sup>* < 1. However, the most notable issue is that the max values reduce as the vegetation patch disappears. Comparing the results of a simple normality test to *ErrorCDF*, it can be concluded that a threshold around *ErrorCDF* = 0.2 can be assumed for the lower limit of intermittency. This threshold is also confirmed by the results of the intermittent turbulence and coherent Reynolds shear stresses at a depth of 0.8 < *z/H<sup>v</sup>* < 1. Correspondingly, in case 4, where *ErrorCDF* < 0.2 at all depths, no anisotropy was observed in the turbulent spectra. In contrast, the anisotropic turbulence was limited to the mixing layer in cases 1, 2 and 3. In the smaller *z/Hv*, the no intermittency is observed for all cases and the effect of bed roughness is the determining factor in the near bed region. Generally, in the near bed region of gravel bed channels, the effect of bed roughness exceeds the patch-produced vortexes [53–55].

**Table 2.** Results of the normality test of velocity time series of case 1 at a distance of *x/Hv*=8. Note that the unmatched values are highlighted.


According to the results of this research, the intermittency and anisotropy of the turbulence are associated with the non-Gaussian distribution of instantaneous velocity. This relationship was also identified by researchers, and is known as a fundamental characteristic of the anisotropic intermittent turbulence [49–51]. However, in the present research, this characteristic was detected in the mixing layer beyond vegetation patches with a relatively high blockage ratio (*Dv/D<sup>c</sup>* > 0.5).

**Figure 9.** Values of *ErrorCDF* for different depths behind the decaying patch at a distance of *x/Hv* = 8. **Figure 9.** Values of *ErrorCDF* for different depths behind the decaying patch at a distance of *x/H<sup>v</sup>* = 8.

According to the results of this research, the intermittency and anisotropy of the turbulence are associated with the non-Gaussian distribution of instantaneous velocity. This relationship was also identified by researchers, and is known as a fundamental characteristic of the anisotropic intermittent turbulence [49–51]. However, in the present research, this characteristic was detected in the mixing layer beyond vegetation patches with a relatively high blockage ratio (Dv/Dc > 0.5). In addition, the formation of intermittent turbulence in the downstream region of the patch has been investigated in this research. One can see from Table 1, by comparing to a normally distributed data set, that the non-Gaussian data sets are dependent on a greater concentration of data between μ-σ and μ + σ. In addition, the values of the normalized standard deviation are higher in these velocity time series. Consequently, the In addition, the formation of intermittent turbulence in the downstream region of the patch has been investigated in this research. One can see from Table 1, by comparing to a normally distributed data set, that the non-Gaussian data sets are dependent on a greater concentration of data between µ − σ and µ + σ. In addition, the values of the normalized standard deviation are higher in these velocity time series. Consequently, the variation in the probability density functions (PDF) of the instantaneous velocity was considered along the flow direction. The results have been compared with the turbulent spectra of the corresponding points. For the mixing layer of case 1, these results are shown in Figure 10. As can be seen in Figure 10, the intermittent anisotropic turbulence formed at a depth of 8 < *x/H<sup>v</sup>* < 12, which is the same range as the occurrence of the outward quadrant dominancy, leading to strong coherent Reynolds shear stresses and anomalies in streamwise velocity. This event can also be inferred from the flattened shape of PDF, which is associated with high standard deviations.

variation in the probability density functions (PDF) of the instantaneous velocity was considered along the flow direction. The results have been compared with the turbulent spectra of the corresponding points. For the mixing layer of case 1, these results are shown in Figure 10. As can be seen in Figure 10, the intermittent anisotropic turbulence formed at a depth of 8 < *x/Hv* < 12, which is the same range as the occurrence of the outward quadrant dominancy, leading to strong coherent Reynolds shear stresses and anomalies in streamwise velocity. This event can also be inferred from the flattened shape of PDF, which is associated with high standard deviations. Furthermore, the intensity of the intermittency can also be evaluated in the energy domain. Researchers used an analysis method based on the filtering of TKE time series via a certain threshold [56,57]. These thresholds include *MTKE* − *aσTKE* and *M + aσTKE*, where *MTKE* is the mean, "*a*" is a constant multiplayer (commonly 1 and 3) and *σTKE* is the standard deviation of the TKE time series in 2D space. The values that exceed this range are considered as strong occurrences characterizing intermittent events. Although the ratio of filtered samples describes the quantity of intermittent events, it cannot provide enough information about the intensity of the intermittency. It is reported that the TKE of extreme events is a suitable indicator of the intermittency [58]. Since extreme events are associated with coherent occurrences in the mixing layer, the Turbulent Kinetic Energy of these extreme events (TKE<sup>i</sup> ) can be used for this propose. This approach is used in this paper for calculating TKE and TKE<sup>i</sup> in a 3D space. Overall, Figure 11 summarizes the results regarding the formation of intermittent turbulence in the mixing layer behind the vegetation patch by means of demonstrating the filtered time series of the instantaneous velocity and the variation in the normalized standard deviation of TKE and normalized TKE<sup>i</sup> .

**Figure 10.** Probability density function and the turbulent spectra for the mixing layer of Case 1. **Figure 10.** Probability density function and the turbulent spectra for the mixing layer of Case 1.

### Furthermore, the intensity of the intermittency can also be evaluated in the energy *4.5. Transformation of Coherent Structures beyond a Vegetation Patch*

domain. Researchers used an analysis method based on the filtering of *TKE* time series via a certain threshold [56,57]. These thresholds include *MTKE-aσTKE* and *M + aσTKE*, where Considering the results of the previous sections, the transformation of coherent structures beyond a vegetation patch can be summarized as below:


(3) With a patch width ratio of *Dv/D<sup>c</sup>* = 0.5, a new type of coherent shear stresses emerges in the wake zone. These structures grow as the dimensions of the vegetation patch reduce, and they come to cover most area of the wake zone with a patch width ratio of *Dv/D<sup>c</sup>* = 0.33. *Water* **2021**, *13*, x FOR PEER REVIEW 19 of 23

**Figure 11.** Formation of the intermittent fluctuations in the mixing layer beyond the patch (case 1). The filtered velocity time series is shown as intermittent fluctuations above the original velocity time series. Variations in the normalized standard deviation of TKE and normalized intermittent Turbulent Kinetic Energy are shown at the bottom. **Figure 11.** Formation of the intermittent fluctuations in the mixing layer beyond the patch (case 1). The filtered velocity time series is shown as intermittent fluctuations above the original velocity time series. Variations in the normalized standard deviation of TKE and normalized intermittent Turbulent Kinetic Energy are shown at the bottom.

*4.5. Transformation of Coherent Structures beyond a Vegetation Patch*  Considering the results of the previous sections, the transformation of coherent structures beyond a vegetation patch can be summarized as below: 1). In the downstream region of a fully channel-spanning vegetation patch, the coherent structures are observable just behind the patch. These structures originate from the stem-scale vortexes that are formed by the leaking flow passing through the vegetation patch [24,28]. As the patch width ratio of *Dv/Dc* reduces to 0.66, the leaking flow In summary, the flow in the downstream region of a decaying vegetation patch alters as the size of vegetation patch decreases. This change includes not only the intensity of occurrences, but also the nature of structure. This means that the dominant coherent structures are transferred from the intermittent anisotropic turbulent fluctuations in the mixing layer of a fully covered patch into the isotropic turbulence associated with the von Karman vortexes on the horizontal plane of the wake layer in the presence of the smallest patch.

increases. Consequently, the stem-scale coherent structures are observed in a larger area (Figure 5). With the patch width ratio of *Dv/Dc* = 0.5, these structures spread to the near-bed region and create a large area of coherent shear stress, in which a Q1-dominant core is surrounded with a sweep ejection-dominant region. However,

layer behind the patch;

### **5. Conclusions**

In the present study, the changes in turbulent structures in the regions downstream of vegetation patches have been investigated. The model vegetation was selected carefully to simulate aquatic vegetation patches in natural rivers. Velocity profile, TKE, turbulent power spectra and quadrant analysis have been used via various approaches to investigate the features and intensity of the turbulent structures. Three different flow layers were detected in the downstream regions of vegetation patches, including the wake layer, the mixing layer and the shear layer. Overall, the results of the TKE and velocity analysis show that those three flow layers (wake, mixing and shear) occurred in all cases of patch setup. However, the wake and shear layers were considerably affected by the flow passing around the sides of the patches. Consequently, different flow structures and associated length scales can be formed, which are dependent on the upward movement of the center of the vertical rotation beyond the patch. The results of spectral analysis indicate that the strong coherent shear stresses are associated with the peaks of the anisotropic turbulent spectra. In contrast, the weak form of coherent shear stresses occurs in the presence of peak isotropic turbulence. The presence of von Karman vortexes near partially covered patches at a depth of *Hv/H* > 0.55~0.7 shows that the frequency of vortex shedding was around 0.07–0.11 Hz, which is equivalent to an average Strouhal number of about 0.25. The main characteristics of turbulent structures beyond a wide patch are associated with the highly intermittent anisotropic turbulent events in the mixing layer, and appear at a depth of *z/H<sup>v</sup>* = 0.7~1.1 and distance of *x/H<sup>v</sup>* = 8~12. The outward interactions are the representative quadrant class of coherent Reynolds shear stresses in the regions downstream of vegetation patches particularly the strong coherent occurrences. The intensity and extent of these structures decrease as the size of the patch reduces. In addition, the intermittency and anisotropy of the turbulence are associated with the non-Gaussian distribution of the instantaneous velocity detected in the mixing layer beyond patches with a relatively high patch width ratio of *Dv/D<sup>c</sup>* > 0.5. Finally, when the size of the patch reduces, von Karman vortexes appear in the wake layer, and form the dominant flow structures in the downstream region of the vegetation patch. It should be noted that the results of this research are limited by our use of a certain type of vegetation, and so it is strongly suggested that other types of the vegetation patches often found in natural rivers be considered in future research.

**Author Contributions:** M.K.; laboratory works, methodology, software and writing—original draft preparation, H.A.; laboratory supervisory, methodology and validation, J.S.; methodology and writing reviewing. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** This study did not involve humans or animals.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are available upon request.

**Acknowledgments:** Special thanks to G. Goodrati Amiri for his contribution in providing the laboratory devices and Esmaeel Dodangeh for his contribution in laboratory setup.

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

### **References**


**Masoud Kazem <sup>1</sup> , Hossein Afzalimehr 1,\* and Jueyi Sui <sup>2</sup>**


**Abstract:** By using model vegetation (e.g., synthetic bars), vortex structures in a channel with vegetation patches have been studied. It has been reported that vortex structures, including both the vertical and horizontal vortexes, may be produced in the wake in the channel bed with a finite-width vegetation patch. In the present experimental study, both velocity and TKE have been measured (via Acoustic Doppler Velocimeter—ADV) to study the formation of vortexes behind four vegetation patches in the channel bed. These vegetation patches have different dimensions, from the channel-bed fully covered patch to small-sized patches. Model vegetation used in this research is closely similar to vegetation in natural rivers with a gravel bed. The results show that, for a channel with a small patch (Lv/Dc = 0.44 and Dv/Dc = 0.33; where Lv and Dv are the length and width of patch and Dc is the channel width, respectively), both the flow passing through the patch and side flow around the patch have a considerable effect on the formation of flow structures beyond the patch. The results of further analysis via 3D classes of the bursting events show that the von Karman vortex street splits into two parts beyond the vegetation patch as the strong part near the surface and the weak part near the bed; while the middle part of the flow is completely occupied by the vertical vortex formed at a distance of 0.8–1 H<sup>v</sup> beyond the vegetation patch, and thus, the horizontal vortexes cannot be detected in this region. The octant analysis is conducted for the coherent shear stress analysis that confirms the results of this experimental study.

**Keywords:** coherent flow structures; vegetation patch; von Karman vortex; octant analysis
