**About the Special Issue Editors**

**Jian Guo Zhou** graduated from Wuhan University with a BSc in River Mechanics and Engineering and subsequently completed his MSc in Fluvial Mechanics at Tsinghua University. He received his PhD in Fluid Mechanics from the University of Leeds. He specialises in formulating mathematical models and developing numerical methods for flow problems in fluids and water engineering, such as coastal, estuary, hydraulic, river and environmental engineering.

**Jianmin Zhang** graduated from Shihezi University with a BSc, and subsequently completed his MSc at Xinjiang Agricultural University. In 2000, he received his PhD in Hydraulics and River Dynamics from Sichuan University. His core research focuses on engineering hydraulics, flood discharge, energy dissipation, hydraulic experiments, two-phase flows and numerical simulations of turbulence.

**Yong Peng** graduated from Northwest A&F University with a BSc in Hydraulic and Hydropower Engineering, and subsequently completed his MSc in Hydraulics and River Dynamics at Sichuan University. He received his PhD in Computational Hydraulics from the University of Liverpool. His main research focuses on flooding, contaminant transport, sediment transport, the lattice Boltzmann method, shallow water flow, two-phase flow, cavitation, and hydraulic experiments.

**Alistair G. L. Borthwick** has more than 40 years' experience in civil engineering. He holds a BEng and a PhD from the University of Liverpool, a DSc degree from Oxford University, and an honorary degree from the Budapest University of Technology and Economics. His research interests include environmental fluid mechanics, coastal and ocean engineering, and marine renewable energy. In 2019, he was awarded the Gold Medal of the UK Institution of Civil Engineers.

## **Preface to "Advances in Hydraulics and Hydroinformatics"**

Great progress has been made in the research on hydraulics, hydrodynamics and hydroinformatics over the past few decades. This includes theoretical, experimental and numerical studies, leading to a new understanding and knowledge of water-related problems, covering a wide range of topics and applications such as hydrology, water quality, river and channel flows. For example, coherent vortex structures in a backward-facing step flow and secondary flow in an open channel bend are measured using PIV; the flow through a Y-shaped confluence channel partially covered with rigid vegetation on its inner bank is measured by ADV. Meanwhile, the velocity decay of an offset jet in a narrow and deep pool, and the skimming flow over a pooled stepped spillway, are studied numerically. In addition, mechanisms for sediment transport due to riverbank failure and a method of predicting maximum scour depth are investigated, and cavitation bubble collapse and its impact on a wall are studied. In order to accelerate knowledge transfer in resolving engineering problems, this Special Issue reports both on the state-of-the-art of the aforementioned research topics, and on advances in sediment transport dynamics, two-phase flows, flow-induced vibration, and hydropower station hydraulics. The goals of this Special Issue are to improve our understanding of the foregoing flow problems and to stimulate future research in these areas, aiming for an improved quality of life.

> **Jian Guo Zhou, Jianmin Zhang, Yong Peng, Alistair G. L. Borthwick** *Special Issue Editors*

### *Article* **Influence of Flushing Velocity and Flushing Frequency on the Service Life of Labyrinth-Channel Emitters**

### **Zhangyan Li 1, Liming Yu 1,\*, Na Li 1, Liuhong Chang <sup>2</sup> and Ningbo Cui <sup>3</sup>**


Received: 19 October 2018; Accepted: 8 November 2018; Published: 12 November 2018

**Abstract:** Dripline flushing is an effective way to relieve emitter clogging and extend the longevity of drip irrigation systems. This laboratory study was conducted at Kunming University of Science and Technology to evaluate the effect of three targeted flushing velocities (0.3, 0.6, and 0.9 m/s) and four flushing frequencies (no flushing, flushing daily, and flushing every three or five days) on the emitter's service life and the particle size distribution of the sediment discharged from emitters and trapped in an emitter channel. The gradation of particle size was analyzed by a laser particle size analyzer. The experiment results suggested that flushing velocity and flushing frequency had a significant effect on the service life of emitters, and the emitter's service life was extended by 30.40% on average under nine different flushing treatments. Flushing can effectively reduce the accumulation of sediments in the dripline and decrease the probability of coarse particles flowing into emitters and fine particles aggregating and cementing in the labyrinth channel, thus relieving the emitter clogging. Therefore, dripline flushing can effectively slow down clogging in muddy water drip irrigation system. The recommended flushing velocity should be set at 0.6 m/s, and the flushing intervals should be shortened.

**Keywords:** clogging; drip irrigation; flushing; particle size distribution; sediments

### **1. Introduction**

Dripline flushing is a necessary maintenance practice for micro-irrigation systems. It removes particles that are not strained by the micro-irrigation system filters and that accumulate in the driplines [1,2]. Physical clogging caused by solid particles is considered the most common emitter clogging category of emitters [3–5]. For drip irrigation systems in the Yellow River irrigation areas of Ningxia and Inner Mongolia, the average sediment concentration in the water abstracted from the Yellow River reaches 35 kg/m3. A large quantity of sand enters into drip irrigation systems and results in emitter clogging even after deposition and prefiltration treatment measures are taken [6]. Therefore, flushing is required to ensure a long economic service life of drip irrigation systems [7]. On the other hand, numerical simulations on shallow water flows and river flushing have been carried out by Peng et al. [8–14].

In order to achieve an optimal flushing effect, drip irrigation systems should be designed properly. Flushing must be done often enough and at an appropriate velocity to dislodge and remove the accumulated sediments [15]. A minimum flushing velocity of 0.3 m/s was recommended by the

American Society of Agricultural and Biological Engineers Engineering Practice, EP-405 [16], but some researchers have advised that a flushing velocity of 0.5–0.6 m/s is necessary when larger particles need to be discharged; for example, it is useful when coarser filters are used in drip irrigation systems [17,18]. In a 30-day field study in which the target flushing velocities were ranged from 0.23 to 0.61 m/s, Puig-Bargués et al. [18] did not find remarkable effects of flushing velocity on the emitter discharge. However, the higher the flushing velocities were set, the more solids were removed from the laterals. Increasing the velocity of flushing may need more costly system designs (e.g., larger supply pipe main, higher pumping requirements, reduced zone sizes) and labor requirements for flushing in irrigation systems should be increased [14].

Different studies explored different flushing frequencies: daily [19], weekly [20], every two weeks [17,21], and monthly [18]. Li et al. [21] claimed that emitter clogging was minimized in cases of biweekly flushing. However, Puig-Bargués et al. [22] reported no obvious differences between flushing frequencies at a velocity of 0.6 m/s, and emitter clogging was primarily affected by the interactions between emitter type, emitter location, and frequency of flushing. Thus, a general agreement on the optimum flushing frequency is lacking at this time, and related studies on flushing frequency and flushing velocity are scarce.

Therefore, the objective of this study was to analyze the influence of three flushing velocities and four flushing frequencies on emitter clogging in drip irrigation systems when using muddy water with full tests. Additionally, particle size distributions of the discharged sediments and residual sediments in the emitter were investigated to determine appropriate flushing schemes to extend the service life of drip irrigation systems.

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

### *2.1. Emitter Characteristics*

The drip tape with non-pressure-compensating emitter manufactured by Dayu Water Conservation Ltd. in Jiuquan City, Gansu Province, China, which is widely used in agricultural irrigation, was applied in the experiments. Its external diameter was 16 mm, and the thickness of the wall was 0.40 mm. Figure 1 presents the structure of the labyrinth-channel emitter. The working pressure was 0.1 MPa. The rated discharge (q) was 1.37 L/h. The flow path had a width (W) of 0.94 mm, length (S) of 37.8 mm, and depth (D) of 0.60 mm. The sectional area (A) was 0.56 mm2. The angle between the bevel edge and baseline was 38.2◦. Space width (L) and height (H) of the tooth were 1.50 and 0.88 mm, respectively. The manufacturing variation of the emitter was 0.028. The flow exponent (x) was 0.49, and the discharge coefficient (k) was 1.12.

The test driplines were cut randomly from the same roll of irrigation drip tapes. The emitter's discharge was measured using pressures of 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, and 0.10 MPa to calculate the discharge coefficient and flow exponent.

**Figure 1.** *Cont.*

**Figure 1.** Emitter used in drip tape. (**a**) Labyrinth-channel emitter; (**b**) parameters of the flow path in the emitters. W, tooth space width; D, depth; α, angle of the bevel edge and baseline; L, tooth space width; h, tooth height.

### *2.2. Experimental Setup*

Figure 2 shows that clean water was stored in a 100-L tank and muddy water was stored in a 150-L tank. Each was equipped with a 1.8-kW submersible pump with a 42-m rated head and 1.8-m3/h rated discharge to provide the working pressure. The pressure adjustment valve and manometer I were installed on the three parallel tubes, which were connected to the submerged pump in the clean water tank, to adjust the different targeted flushing velocities. Manometer II was used to monitor the irrigation pressure. The pressure in manometer I and manometer II was 0.06 and 0.25 MPa, respectively, and their precision was 0.25%. A test platform (1.5 m in width, 4.8 m in length, and 1.0 m in height) was applied to support the laterals and emitters. Ten driplines were arranged on the test platform and all of them were equipped with control valves at the front and back ends. The spacing between each lateral was 0.16 m. Each lateral pipe had 15 emitters, with a spacing of 0.3 m.

**Figure 2.** Experimental setup. 1, muddy water tank; 2, submersible pump; 3, clean water tank; 4, manometer I; 5, manometer II; 6, pressure adjustment valve; 7, polyvinyl chloride (PVC) pipe; 8, control valve; 9, PVC choke plug; 10, lateral; 11, measuring cup.

### *2.3. Measurement of Particle Size Distribution and Water Source*

The distribution of particles size was analyzed by the Malvern laser particle size analyzer 2000 (Malvern Instruments Ltd., Malvern, UK). The measuring range of the particle analyzer was 0.02–2000 μm, and the particles were arranged in order of increasing size. When the accumulation volumes of sediment reached 10%, 50%, and 90%, the largest values of particle size were noted as D10, D50, and D90, respectively.

Tap water was used for lateral flushing in this study. Based on the maximum sediment concentration of irrigation water (0.8 g/L), the sand content in the muddy water was set at 2.5 g/L to accelerate the clogging process and reduce the experimental period. The test sediment was sandy loam soil obtained from Kunming, Yunnan Province of China, which was filtered using screens with 160 meshes after natural air drying. (In this study, 120 mesh screens were not selected for filtration to avoid insufficient test numbers caused by untimely clogging and the deposition of coarse particles.) Five sediment samples were randomly taken from the filtered sand after mixing well and the test sediment particle size distribution was presented based on their average values (see Table 1). The cumulative distribution and differential distribution of the tested sediments are shown in Figure 3.

**Figure 3.** Differential distribution and cumulative distribution of sediments.

**Table 1.** Size distribution of sediment particles.


### *2.4. Test Procedures*

An experimental setup was conducted in a laboratory at Kunming University of Science and Technology, Kunming, Yunnan Province of China. In order to prevent the emitter flow rate and anti-clogging performance from being influenced by temperature variation, the experiments in this study were conducted between 13:00 and 16:00 CST, from 19 April to 7 September in 2017. Repeat experiments were carried out in three periods due to the limited width of the platform. The first repeated test was conducted between 19 April and 2 June, the second between 5 June and 19 July, and the third between 25 July and 7 September. The average temperature of the environment for the three test periods was 21.39 ◦C, 22.27 ◦C, and 22.08 ◦C, respectively.

### 2.4.1. Muddy Water Irrigation

The test conditions were as follows: normal irrigating pressure, 0.1 MPa; duration of irrigation event of each test, 20 min; initial volume of muddy water, 100 L; and total volume for each irrigation event, approximately 68.5 L. The muddy water was well stirred manually to avoid sedimentation. The emitter's discharge was measured every day using a marked 1000-mL measuring cup which was placed exactly beneath all the emitters. After each irrigation event, the weight of each measuring cup containing emitted muddy water was measured by a digital balance. The resolution of the digital balance was 0.01-g.

The emitter discharge was recorded as volume per unit time (L/h). For the emitter discharge less than 75% of the nominal rated flow, it was deemed to be clogged. The tested muddy water was replaced with new muddy water every day.

### 2.4.2. Sampling and Testing

The irrigation test was conducted with 45 events and 10 driplines, and the average discharge of all emitters in each lateral was taken into consideration. Therefore, 10 samples of each irrigation were measured three times. In all the tests, a total of 1350 (i.e., 45 × 10 × 3 = 1350) muddy water samples were evaluated. The distribution of particle size was measured with all emitters in each dripline. Ten samples were assigned to each irrigation event, and the process was repeated three times. A total of 1350 samples were measured.

### 2.4.3. Clean Water Flushing

The present study evaluated the effect of two factors: (1) the flushing frequency, either daily (F1), once per three days (F1/3), or once per five days (F1/5); and (2) the flushing velocity, at either 0.3 m/s (V0.3), 0.6 m/s (V0.6), or 0.9 m/s (V0.9). Table 2 shows 10 treatments consisting of a control non-flushed treatment (T10) and a combination of targeted flushing velocities and flushing frequencies (T1–T9). There was only one dripline being flushed in every experiment, and each flushing event lasted for 5 min.

**Table 2.** Flushing velocity and flushing frequency treatment.


### *2.5. Statistics*

The data were analyzed using SPSS software for Windows, version 22.0 (IBM Corp., Chicago, IL, USA). Based on main effects analysis of variance (ANOVA), the significance of the differences between treatments for different response variables was evaluated. Furthermore, the results were categorized as "descriptive statistics" (such as mean) and Levene's test of equality of error variances. In addition, multiple comparisons were carried out based on Fisher's least significant difference tests because significant differences (*p* < 0.05) were indicated by ANOVA.

### **3. Results**

### *3.1. Variations in Emitter Discharge*

Figure 4 illustrates the value changes of the emitter discharge versus irrigating events. The level of straight line was 75% of the initial discharge, which was used as the criterion for emitter inefficiency or severe emitter clogging. Ten discharges of the emitters treated under diversified conditions showed a descendant trend with increasing irrigation events, suggesting that clogging with different degrees was seen in emitters until their service life ended. In the figure, the bold black line shows changes of the emitter discharge in T10, the control experiment. The emitter discharge in treatment T10 reached 75% of the rated discharge first, and its average discharge at the end of the test was the minimum without any flushing treatments made after irrigation. According to the irrigating events (Table 2), T10 provided only 23.33 irrigation events before emitters were severely plugged. However, for the remaining nine flushing treatments, the average normal irrigating event was 33.52, which indicated that the emitter's average service life in the flushed driplines after irrigation increased by 30.40% in

comparison with those in the non-flushed T10. In the nine flushing treatments, treatment T7 had the longest irrigation event 41.33 and the universally largest average discharge; treatment T3 had the shortest irrigation event 27.33 and the universally lowest average discharge after 28 irrigation events. The service life of treatments T7 and T3 increased by 43.55% and 14.64%, respectively, compared with the non-flushing group T10. Additionally, treatments T7 and T3 were flushed at the frequency of F1 and F1/5 and at the velocity of V0.9 and V0.3, respectively, demonstrating that the emitters had a longer service life under high flushing velocities and high flushing frequencies.

**Figure 4.** Average change in the emitter discharge versus irrigating events.

### *3.2. Variance Analysis of the Two Flushing Factors on the Service Life of Emitters*

Flushing velocity and flushing frequency have significant effects on the service life of emitters, and the interactions were negligible (see Table 3). Figure 5a shows that the service life of emitters increased with the increase in the flushing velocity. The average service life of emitters was 22.80%, 31.60%, and 35.59% higher under V0.3, V0.6, and V0.9 conditions, respectively, compared with the non-flushing conditions. Moreover, the service life of emitters under V0.6 and V0.9 conditions was similar, which meant that the flushing velocities in excess of 0.6 m/s would only help improve the service life to an insubstantial degree. Therefore, 0.6 m/s proved to be the optimal flushing velocity.

**Table 3.** Variance analysis of the effects of flushing velocity and flushing frequency on the service life of emitters.


\* *p* < 0.05; \*\* *p* < 0.01.

**Figure 5.** Average service life and standard errors of emitters at different flushing (**a**) velocities and (**b**) frequencies.

On the contrary, the service life of emitters was significantly reduced as the flushing frequency decreased (Figure 3b). The average service life of emitters was 38.43%, 30.71%, and 19.55% higher under F1, F1/3, and F1/5 conditions, respectively, compared with the non-flushing treatment, indicating that a decrease in the flushing frequency effectively extended the service life of emitters.

### *3.3. Particle Size of Discharged Sediments*

Table 4 shows that the effect of flushing velocity on D90 was significant, indicating that coarse particles passing through the labyrinth channel were affected by the flushing velocity, whereas fine particles were not. The flushing frequency had no effect on the particle size of discharged sediments. Figure 6 shows that D10, D50, and D90 mean values under non-flushing measures were 3.00, 7.16, and 12.65 μm, respectively, which were significantly higher than the different processing levels of flushing velocity and flushing frequency. This can be attributed to the fact that residual sediments in the lateral were discharged by flushing, resulting in a reduced concentration of sediments at the inlet of the labyrinth channel. Figures 7 and 8 show large amounts of residual sediments in the non-flushing lateral. The sediments mixed together, and more coarse particles flew into the labyrinth channel with the turbulence of the flow. Therefore, flushing effectively reduced the concentration of residual sediments in the lateral and decreased the probability of large particles entering into the labyrinth channel. Meanwhile, Figure 6b shows that the degree of D90 increased slightly with the decrease in the flushing frequency for the same reason.

**Table 4.** Variance analysis of the effects of flushing velocity and flushing frequency on the particle size of discharged sediments.


<sup>\*</sup> *p* < 0.05; \*\* *p* < 0.01.

**Figure 6.** Average and standard errors of D10, D50, and D90 at different flushing (**a**) velocities and (**b**) frequencies.

### *3.4. Residual Sediments in Drip Tape and Emitters*

Emitters and pipe sampled from the driplines were cut open to provide visual proof of clogging when experiments were finished (Figures 7 and 8). A few sediments were observed on the inner walls of the lateral with flushing treatment, whereas substantial sediment accumulations were observed in the lateral without flushing (Figure 5). Hence, flushing can effectively prevent particles from accumulating in the drip tape and emitters. Figure 8 shows that the sediment particles were seldom deposited in the labyrinth channel and the outlet of emitters due to flow turbulence with flushing, whereas clogging particles were found along the labyrinth channel of emitters without flushing, especially at the outlet. Additionally, other impurities such as plant root debris, which were not or seldom observed in the labyrinth channel of emitters with flushing treatment, were widely distributed in the labyrinth channel of emitters without flushing. In summary, flushing could effectively remove residual sediments and other impurities in the labyrinth channel, thus enhancing the anti-clogging performance of emitters.

**Figure 7.** Residual sediments in the drip tape after tests. (**a**) Drip tape with flushing treatment; (**b**) drip tape without flushing treatment.

**Figure 8.** Residual sediments in the labyrinth channel of emitters (yellow parts indicate that the clogging substances are other impurities, not sediments). (**a**) Emitters with flushing treatment; (**b**) emitters without flushing treatment.

Table 5 shows the variance analysis of clogging substances in emitters. The flushing frequency and flushing velocity had significant effects on D10, D50, and D90. Figure 9 shows that D50 and D90 were the highest under non-flushing conditions, being 30% and 52% higher than the averages under the corresponding flushing conditions, respectively, indicating that flushing prevented large particles from entering into the labyrinth channel. The figure also shows that D50 and D90 decreased as the flushing velocity increased. For instance, D10, D50, and D90 under V0.9 conditions were the minimum, being, respectively, 25.76%, 38.55%, and 47.14% lower than those under non-flushing conditions. In summary, a higher flushing velocity contributed to the discharge of coarse particles in emitters, thereby lowering risks of emitter clogging.

**Table 5.** Variance analysis of the effects of flushing velocity and frequency on the particle size of residual sediments in emitters.

**Figure 9.** Average and standard error of D10, D50, and D90 at different flushing (**a**) velocities and (**b**) frequencies.

(**a**) (**b**)

Figure 9 further shows that D50 and D90 increased as the flushing frequency decreased and was maximized under non-flushing conditions. D90 under F1, F1/3, and F1/5 conditions was 48.50%, 36.63%, and 17.66%, respectively, significantly lower than that under non-flushing conditions. D50 under F1 and F1/3 was 37.60% and 25.85%, respectively, significantly lower than that under non-flushing conditions. However, D10 differed slightly between the F1, F1/3, and F1/5 conditions. A low flushing frequency caused the accumulation of sediments in the lateral and increased the probability of large particles being trapped in the labyrinth channel, resulting in increased particle sizes of residual sediments in the channel. It also favored the agglomeration of sediments, resulting in a higher degree of D50 and D90. Combining these two effects, the flushing frequency had a significant effect on the size distributions of residual sediments in the labyrinth channel.

### **4. Discussion**

### *4.1. Influence of Flushing Treatment on Emitter Clogging*

For drip irrigation systems, appropriate flushing can effectively prevent emitter clogging by hindering the agglomeration of sediment particles or their adhesion to organic residuals to generate large particles [23,24]. This study demonstrated that lateral flushing using clean water after the irrigation of muddy water could enhance the service life of emitters by 14.64–43.55% compared with non-flushing cases, and both the flushing velocity and flushing frequency significantly affected the service life of emitters. However, Puig-Bargués et al. [18] reported that the flushing frequency and flushing velocity had little effect on the resulting emitter discharges in subsurface drip irrigation systems measured at the end of the study. The results were different from the findings of the present study mainly due to different irrigation systems and flushing operations.

Recent studies found that the service life of emitters increased with the increase in the flushing velocity. For instance, the service life of emitters flushed at the velocities of 0.6 and 0.9 m/s was extended by more than 30% compared with that without flushing. This paper held that the flushing velocity of 0.6 m/s would appear to be adequate. In a laboratory study simulating a dripline with a transparent PVC pipe, Puig-Bargués and Lamm [25] suggested that the minimum flushing velocity of 0.3 m/s appeared sufficient for most micro-irrigation systems working under typical conditions. This was not consistent with the findings in this paper, which could be attributed to the fact that the researchers did not fully consider the complexities of the flow regime that might occur in the emitter. Zhang et al. [26] compared the influence of constant pressure and pulse pressure on the anti-clogging performance of the labyrinth emitter. According to their findings, the pulse pressure could alleviate clogging and reduce discharge of the emitter. Some studies [27,28] demonstrated that the emitter's anti-clogging performance could be enhanced by increasing the flushing pressure. This was consistent with the conclusion of this study because the pressure of flushing was proportional to the velocity of flushing.

In this study, the service life of emitters degraded significantly in the order of F1 > F1/3 > F1/5 > without flushing. Li et al. [21] conducted an in situ surface drip irrigation experiment under three conditions of dripline flushing frequency, in which the recycled water from a sewage disposal plant was used. They reported that lateral flushing could significantly relieve the emitter clogging of the irrigation systems using reclaimed water. Elberry et al. [23] reported that clogging at a high flushing frequency was significantly lower than that at a low flushing frequency. The aforementioned study results are consistent with the findings of this study. Nevertheless, Feng et al. [29] found that emitter clogging was not prevented with five different frequencies of flushing in the drip irrigation experiment using salty groundwater, which is primarily affected by flushing operations and water quality.

### *4.2. Influence of Flushing Treatment on Particle Size Distributions*

Flushing can effectively enhance the anti-clogging performance of emitters by removing residual sediments from the dripline and maintaining the cleanliness and the passability of the labyrinth channel [29,30]. The results of this study indicated that the passability of large sediment particles in the labyrinth channel was the dominant factor for emitter clogging. In this study, the corresponding pressure in emitters was 0.004, 0.008, and 0.016 MPa, which was much less than 0.1 MPa, the normal irrigation pressure under the flushing velocity of 0.3, 0.6, and 0.9 m/s, respectively. Nevertheless, this study demonstrated that these relatively low pressures still had a significant effect on the particle size of residual sediments in the labyrinth channel. Kou et al. [31] reported that a disordered system of granular materials, such as powders, sand, and foams, formed stable structures when unperturbed, but that they "relax" in the presence of external influences such as shear or tapping, becoming fluid in nature. Similarly, the slightest perturbation of a particle system due to flushing destabilized some of the particles trapped in the labyrinth channel of emitters, increasing the flowability of particles in the channel. Hence, the passability of the labyrinth channel was enhanced, thus extending the service life of emitters. In this study, D90 of discharged sediments was significantly affected by the flushing velocity. Specifically, D90 under V0.9 was significantly larger than that under V0.3 and V0.6 and without flushing, indicating that the perturbation of the particles increased with the increase in the flushing velocity and resulted in the increased trafficability of coarse particles.

A certain concentration of particles is indispensable for the follow-up collision and flocculation of particles. As particle concentration in the influent water increased, the local concentration distribution of particles in the labyrinth channel was improved [32–34]. In this study, the particle size of discharged sediments increased with the increase in the flushing frequency and reached its maximum under non-flushing conditions. As the flushing frequency decreased, fine particles were continuously discharged from the emitters, resulting in an increase in the relative contents of coarse particles and concentration of sediments, thereby improving the probability of collision, sedimentation, and accumulation of large sediment particles. This was why D50 and D90 of residual sediments in emitters increased with the increase in the frequency of flushing.

### *4.3. Influence of Particle Size on Emitter Clogging*

The clogging mechanism in the emitter is dependent on the particle size of sediments. Coarse particles readily underwent collision and deposition in the vortices of the water flow due to bigger sizes and the larger drag force of the particles. Therefore, the particles could not escape from the vortices. On the other hand, the major forces causing the fine particles to clog came from their surface tension and adhesion to surrounding substances. In this study, D10 and D50 of discharged sediments were not affected by flushing velocity or flushing frequency, indicating that fine particles had better flow characteristics due to their lower drag force in the flow. In other words, this study demonstrated that coarse particles tended to be trapped in the labyrinth channel, which was consistent with the findings of previous studies [35–38].

Wu et al. [39] pointed out that it was difficult to cause clogging when the particle size was less than 20 mm. Particles <20 mm took up 50% of the muddy water (see Table 1), whereas D90 of discharged sediments was significantly lower than 20 μm (Table 6). Moreover, D100 of discharged sediments was only 27.18 μm. This indicated that particles larger than 20 μm tended to stay in the labyrinth channel and that particles smaller than 20 μm could easily flow out with the water. Table 6 shows that D50 of residual sediments in emitters was 24.38–39.67 μm, which was 22–99% higher than that in the initial water. In addition, D90 of residual sediments in emitters was 44.74–86.89 μm, which was 13–120% higher than that in the initial irrigation water resource.

Generally, the agglomeration and flocculation of fine particles occurred during the flow, and the upper-limit particle size was 10 μm [40]. Table 6 shows that the percentage of particle sizes below 10 μm was 10% in the residual sediments in the labyrinth channel but 90% in the discharged sediments. Hence, the emitter plugging in this study was caused mainly by large residual sediment particles in the emitter's labyrinth channel, instead of the agglomeration and flocculation of small sediment particles. This could be attributed to high flushing velocity and high flushing frequency.


**Table 6.** Average particle size of residual sediments in emitters and discharged sediments (μm).

### **5. Conclusions**

The nine flushing measures led to an enhancement of 30.40% in the service life of emitters on average, and both flushing velocity and flushing frequency had significant effects on the service life of emitters. Although the life of emitters flushed at 0.9 m/s was longer than that of emitters flushed at 0.6 m/s, the numerical differences were small (<4%). Therefore, the flushing velocities setting at around 0.6 m/s appeared to be sufficient. Increasing the flushing frequency would be an economical way to achieve the optimal flushing effect without providing a higher flushing velocity, which might increase the system cost. Moreover, the inner walls of the lateral were clean in the presence of flushing and only local clogging was observed in the labyrinth channel. Severe sedimentation was observed in the lateral in the absence of flushing, and the clogging by sediments was observed along the labyrinth channel. Furthermore, flushing decreased the accumulation of sediments in the lateral, which prevented large sediment particles from entering into the labyrinth channel and lowered the probability for small sediment particles to agglomerate and flocculate into large particles, thus reducing the risks of emitter clogging.

**Author Contributions:** Software, Z.L.; Methodology, N.L.; Writing—original draft preparation, Z.L.; Writing review and editing, L.Y.; Funding acquisition, L.Y., L.C., and N.C.

**Funding:** This research was funded by the National Natural Science Foundation of China (NSFC) (Nos. 51769009, 51379024, 51809022) and the Funding Program for University Engineering Research Center of the Yunnan Province in China.

**Acknowledgments:** The authors would like to express their appreciation for the financial support received from the National Natural Science Foundation of China (NSFC) (Nos. 51769009, 51379024, 51809022) and the Funding Program for University Engineering Research Center of the Yunnan Province in China.

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

### **References**


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

### *Article* **Effect of Particle Size and Shape on Separation in a Hydrocyclone**

### **Zhaojia Tang 1, Liming Yu 1,\*, Fenghua Wang 1, Na Li 1, Liuhong Chang <sup>2</sup> and Ningbo Cui <sup>3</sup>**


Received: 24 October 2018; Accepted: 18 December 2018; Published: 21 December 2018

**Abstract:** Given the complex separation mechanisms of the particulate mixture in a hydrocyclone and the uncertain effects of particle size and shape on separation, this study explored the influence of the maximum projected area of particles on the separation effect as well as single and mixed separations based on CFD–DEM (Computational Fluid Dynamics and Discrete Element Method) coupling and experimental test methods. The results showed that spherical particles flowed out more easily from the downstream as their sizes increased. Furthermore, with the enlargement of maximum projected area, the running space of the particles with the same volume got closer to the upward flow and particles tended to be separated from the upstream. The axial velocity of the combined separation of 60 μm particles and 120 μm particles increased by 25.74% compared with that of single separation of 60 μm particles near the transition section from a cylinder to a cone. The concentration of 60 μm particles near the running space of 120 μm particles increased by 20.73% and those separated from the downstream increased by 4.1%. This study showed the influence of particle size and maximum projected area on the separation effect and the separation mechanism of mixed sand particles in a hydrocyclone, thereby providing a theoretical basis for later studies on the effect of particle size and shape on sedimentation under the cyclone action in a hydrocyclone.

**Keywords:** CFD–DEM coupling; hydrocyclone; particle shape; particle size; water and sediment separation

### **1. Introduction**

A hydrocyclone is widely used for various separation purposes in different fields, such as mining, chemical industry, and agriculture. In the engineering field, it is often adopted to separate and classify the sand in water. The specific working process is as follows. As the water and sand mixture enters from the inlet, the water flows out from the upstream and downstream. The sand particles with strong water-following ability are carried out from the upstream and those with poor water-following ability are carried out from the downstream under the internal swirling flow.

In recent years, the separation mechanism of hydrocyclone has been mainly studied from the following aspects:

(1) The separation effect of various solid particle properties and different external conditions: The equation of particle motion in a hydrocyclone is established and the formula for calculating the optimum separation medium density is deduced according to the force analysis on different particles in the separation [1]. The calculation method and equation are put forward based on the relationship between particle size and fluid drag force [2]. The separation effect of submicron particles in a microclone at different temperatures and water pressures [3], the effect of internal flow field on the separation of different-sized particles [4], and the effect of inlet velocity and particle shape on the separation effect [5] are also studied.

(2) The effect of a hydrocyclone structure on separation: The optimal particle separation effect is determined by the orthogonal analysis on the optimal insertion depth, wall thickness, and optimal inlet flow [6]. The influence of the inlet structure in the performance of a hydrocyclone is explored [7]. Three kinds of hydrocyclones are designed to improve the particle separation efficiency, and the optimum separation structure is determined by studying the pressure and velocity distribution inside a hydrocyclone [8]. The effect of cone angle structure on the classification of fine particles [9] and of the inlet angle, cylinder height, and inner cone on separation [10,11] are discussed. In addition, studies are conducted on the separation results with hydrocyclones of diverse structures and different working conditions. Factors such as feed flow rate, size of the downstream, isolation tank, suspension concentration, viscosity, and pressure drop are included [12].

The separation effect of solid particles has been investigated by observing changes in the external conditions, optimization of hydrocyclone structure, and so on. Only a few basic studies have been conducted on the shape and size of sand particles. Thus, more efforts are called for on the study of the interaction between particles and the complex motion of different-shaped particles in a hydrocyclone. Due to the differences in particle sizes and shapes in nature, the separation becomes more complex, with unpredictable separation effects. In this study, the CFD–DEM coupling method and experimental testing method were employed to study the effect of maximum projected area on the separation performance with the same volume of sand. The mechanism underlying the separation of water and sand in a hydrocyclone was proposed by comparing the separation of single 60 μm spherical sand particles and 60 and 120 μm mixed particles.

### **2. Experimental Modes and Methods**

### *2.1. Hydrocyclone Model*

The parameters of the hydrocyclone used in this study were as follows. The diameters of the inlet, upstream, and downstream are 2, 3.4, and 2.2 mm, respectively. The cylinder section was 40 mm in height and 10 mm in diameter. The cone angle was 4 degrees. The shape and mesh conditions are shown in Figure 1. Based on the separation characteristics of hydrocyclones, it is generally believed that the sand particles moving against the wall where the cylinder and cone connect can be separated from the downstream. Therefore, the D-D section was set on the model to divide the area between the radial center and the wall into five regions. Figure 1a shows the model segmentation and meshing, and Figure 1b provides a magnified view of the D-D section 1:5. The three-dimensional structure was drawn using the UG10.0 (Unigraphics NX10.0) software, the meshing was done using the GAMBIT 2.4.6 software, and the numerical simulation was based on the ANSYS 16.0 and DEM 2.6.1 coupling method.

**Figure 1.** The axial length of *Z* axis and cross section of the cylinder of the hydrocyclone.

### *2.2. Numerical Simulation*

### 2.2.1. Particle Characteristics

In this study, a single separation simulation of five spherical sand particles with diameters of 60, 70, 90, 106, and 120 μm was performed. A single separation simulation was conducted on three plate-like particles with the same volume as spherical sand particles of 60 μm (in Figure 2, the maximum projected area is S), and their maximum projected area was 3.11 S (type A), 2.24 S (type B), and 1.36S (type C), respectively. A mixed separation simulation was conducted on spherical sand particles with diameters of 60 and 120 μm. Table 1 shows the deformation amount, thickness, maximum projected area, and equivalent diameters of particles are shown in Table 1. The type C, B, and A particles share the maximum projected area with particles of 70 μm, 90 μm, and 106 μm, respectively.

**Figure 2.** Particle shapes in numerical simulation.

**Table 1.** Parameters of particles in the numerical simulation.


In simulation experiments, the recorded residence time of sand particles in the cyclone equaled the period from their entering the inlet to their leaving the downstream or upstream. The data of particle velocity, position, and time were exported from the software DEM and edited by the MATLAB 2017b (Matrix Laboratory) postprocessor program. The values between the third and fourth seconds were collected as the analysis data.

### 2.2.2. Mathematical Model and Simulation Method

Flow in the hydrocyclone is considered to be an incompressible viscous fluid. The effect of gravity and the roughness of the hydrocyclone walls were taken into account, though surface tension was neglected. In the present study, the inlet velocity of the hydrocyclone was 2.00 m s<sup>−</sup>1, The continuous medium flow is calculated from the continuity and the Navier–Stokes equations based on the local mean variables over a computational cell, which are given by:

$$
\frac{
\partial \varepsilon}{
\partial t} + \nabla \cdot (\varepsilon \mathbf{u}) = 0 \tag{1}
$$

$$\frac{\partial(\rho\_{\text{f}}\varepsilon\mathbf{u})}{\partial\mathbf{t}} + \nabla \cdot (\rho\_{\text{f}}\varepsilon\mathbf{u}\mathbf{u}) = -\nabla\mathbf{P} - \mathbf{F}\_{\text{P}} + \nabla \cdot (\varepsilon\mathbf{r}) + \rho\_{\text{f}}\varepsilon\mathbf{g} \tag{2}$$

Fp−<sup>f</sup> is Interaction forces between fluid and solid phases, equal to kc ∑ i=1 fP−f,i/ΔVc, N/m3; the flow solved in Equations (1) and (2) represents the mixture flow of medium and sand, and was obtained by use of the models in a commercial ANSYS software package, i.e., Fluent. The details of the medium flow calculation and its validation can be found elsewhere [13,14]. In this work, we only give an overall description of the ANSYS model for the hydrocyclone.

The governing equations that are described above were discretized by the control volume numerical technique, and then the SIMPLE pressure–velocity coupling technique, with a second-order upwind scheme for the convection terms that were employed to solve the discretized equations over the computational domain.

The maximum particle volume fraction for the sand used in this study was less than 1%, which means the mixture of water and sand belonged to dilute phase flow. The Lagrangian coupling method was employed and Table 2 summarizes the parameter settings for the sand [15], Table 3 summarizes the boundary condition in the model. Ansys16.0 and DEM2.6.1 were employed during this study. The discrete approach was used to simulate sand movement, collisions among sand particles and between the sand particle and the hydrocyclone wall, and the effects of sand movement on the surrounding continuous phase, energy and momentum exchange. Collisions among sand particles and between sand particles and the wall did not lead to significant plastic deformation, which in consequence were attributed to hard particle contact, which is a wet grain contact model. The 'Hertz-Mindlin (no slip) built-in' model was utilized in this work [16].





In the simulation, the contributions of viscous drag force, pressure gradient force and gravity were considered. Other additional forces such as virtual mass force and Saffman force were not considered, as they are an order of magnitude smaller compared with the foregoing [17].

The CFD-DEM coupling process is described as follows: the continuous phase is resolved by ANSYS to acquire fluid drag force with sand, which was transformed from flow field information through the drag force model. The stress state of sands was measured by DEM to obtain new information such as the position and velocity of sands and to the flow field. ANSYS was used to model the flow and the most representative stress condition of sediments. The two approaches were coupled with a compound model that included the particle–fluid interaction force [16]. The relevant equations can be seen in Table 4.


**Table 4.** Equations and specifications.

Particle–fluid interaction, force (fp−f, i):

ε = 1 −

$$\begin{aligned} \text{f}\_{\text{D,i}} &= \left( 0.63 + \frac{4.8}{\text{Re}\_{\text{p},i}^{0.5}} \right)^2 \frac{\rho\_i |\mathbf{u}\_i - \mathbf{v}\_i| (|\mathbf{u}\_i - \mathbf{v}\_i|)}{2} \frac{\pi \mathbf{d}\_i^2}{4} \varepsilon\_i^{-\beta} \\ \text{Re}\_{\text{p},i} &= \frac{\mathbf{d}\_i \rho\_i \varepsilon\_i |\mathbf{u}\_i - \mathbf{v}\_i|}{\mu\_i}, \ \beta = 3.7 - 0.65 \exp\left[ -\frac{\left( 1.5 - \log \text{Re}\_{\text{p},i} \right)^2}{2} \right] \\ &\stackrel{\text{h}\_C^{\text{C}}}{\sum} \mathbf{v}\_i \end{aligned} \tag{3}$$

$$\frac{\partial \mathbf{V}}{\partial \mathbf{V}\_{\mathbf{C}}}$$

fpg,i = Vp,i∇P (4)

The free settling velocity (ut)

$$\mathbf{u}\_{\rm t} = \frac{\mathbf{g} \mathbf{d}^2}{18\nu} \left( \frac{\rho\_{\rm p} - \rho}{\rho} \right) \tag{5}$$

Here, the drag coefficient (C) is

$$\mathbf{C} = \frac{(\rho\_\mathbf{P} - \rho)\mathbf{V}\mathbf{g}}{\mathbf{A}\rho\left(\frac{\mathbf{u}^2}{2}\right)}\tag{6}$$

And the particle Reynolds number (Re) is

$$\text{Re} = \frac{\text{D}\_{\text{1p}} \text{up}}{\mu} \tag{7}$$

The drag coefficient varies depending on the inverse of the particle Reynolds number, thus the approximated drag coefficient was calculated as a function of the particle Reynolds number. The approximated drag coefficient Ca is expressed as

$$\mathbf{C\_{a}} = \mathbf{C\_{0}} + \mathbf{C\_{1}} \left(\frac{1}{\text{Re}}\right) \tag{8}$$

and C0 and C1 were determined from the correlation of the particle Reynolds numbers and measured drag coefficients using the ratio of the particle diameter D2 to thickness T as the particle shape factor

$$\mathbf{C}\_{0} = 0.4555 \text{Ln} \left( \frac{\text{D}\_{2}}{\text{T}} \right) + 0.4687 \tag{9}$$

$$\mathbf{C}\_{1} = 19.285\tag{10}$$

Finally, the approximated drag coefficient was expressed as [18]:

$$\mathbf{C\_{a}} = 19.285 \left( \frac{1}{\mathrm{Re}} \right) + 0.4555 \mathrm{Ln} \left( \frac{\mathrm{D\_{2}}}{\mathrm{T}} \right) + 0.4687 \tag{11}$$

Pressure gradient force [19]: Due to the uneven distribution of hydraulic pressure on the surfaces of particles, the force is given by:

$$\mathbf{F\_{P}} = -\frac{\nabla \mathbf{p}}{\rho\_{\rm P}}\tag{12}$$

The total drag force of fluid on sand particles (FD)

$$\mathbf{F\_D = \xi A\_p \frac{\rho \mu^2}{2}} \tag{13}$$

When a particle accelerates in the viscous fluid, it will cause the surrounding fluid to accelerate, too. Due to the inertia effect, the fluid has a reaction force to the particle, which makes the inertia of the particle virtually increase. The related force is given by [2]:

$$\mathbf{F}\_{\rm A} = \frac{1}{2} \frac{\rho}{\rho\_{\rm p}} \left( \frac{\rm du}{\rm dt} - \frac{\rm du\_{\rm p}}{\rm dt} \right) \tag{14}$$

When a particle moves in a flow field with velocity gradient, it will experience a lateral force from the low flow rate side to the high flow rate side, given by:

$$\mathbf{F\_S} = \frac{9.69}{\pi \mathbf{d\_{p}} \rho\_{\mathbf{P}}} (\mu \rho |\mathbf{y}|)^{\frac{1}{2}} \mathbf{C\_{LS}} (\mathbf{u} - \mathbf{u\_{P}}) \tag{15}$$

when a particle rotates at an angular velocity of ω, a lateral force perpendicular to the relative velocity and the rotation axis of the particle is generated from the upstream side to the downstream side, given by:

$$\mathbf{F\_M} = \frac{3}{4} \frac{\rho}{\rho\_\mathbf{p}} \frac{\left(\mathbf{u} - \mathbf{u\_P}\right)^2}{\mathbf{d\_P}} \mathbf{C\_{LM}} \frac{\omega\_\mathbf{r} \times \mathbf{u\_r}}{|\omega\_\mathbf{r}||\mathbf{u\_r}|}, \ \omega\_\mathbf{r} = \omega - \Omega \tag{16}$$

When a particle moves in a viscous fluid, the boundary layer on the surface of the particle will carry a certain amount of fluid. Due to the inertia of the fluid, the fluid cannot immediately follow the acceleration or deceleration of the particle, so this boundary layer is unstable, and the particle will be affected by a time-dependent force, given by:

$$\mathbf{F\_B = \frac{9}{\mathbf{d\_P \rho\_P}} \sqrt{\frac{\rho \mu}{\pi}} \int\_{t\_0}^{t} \frac{(\mathbf{d}u/d\tau) - (\mathbf{d\_P}/d\tau)}{\sqrt{t-\tau}} d\tau \tag{17}$$

In the Lagrangian coordinate system, the motions of particles are predicted based on the forces acting on them. In our studied system, the governing equation of particles can be given by

$$\frac{d\mathbf{u\_{P}}}{dt} = \left(1 - \frac{\rho}{\rho\_{\text{P}}}\right)\mathbf{g} + \mathbf{F\_{D}} + \mathbf{F\_{P}} + \mathbf{F\_{A}} + \mathbf{F\_{S}} + \mathbf{F\_{M}} + \mathbf{F\_{B}}\tag{18}$$

The formula for the separation effect of the hydrocyclone:

$$\begin{aligned} \text{P\_{downstream}} &= \frac{\text{N}\_{\text{downstream}}}{\text{N}\_{\text{downstream}} + \text{N}\_{\text{upstream}}} \% \text{o} \\ \text{P\_{upstream}} &= \frac{\text{N}\_{\text{upstream}}}{\text{N}\_{\text{downstream}} + \text{N}\_{\text{upstream}}} \% \text{o} \end{aligned} \tag{19}$$

### *2.3. Separation Experiments*

Figure 3 provides the flow chart of separation experiment in the hydrocyclone. The water pressure was controlled at 20 kPa by adjusting the feed valve (consistent with the pressure set in numerical simulation).

**Figure 3.** Schematic diagram of the hydrocyclone test.

Sand particles with D50 of 120 μm were screened using a 115 to 125 mesh screen to substitute the 120 μm particles in numerical simulation, and particles with D50 of 60 μm were screened using a 230 to 250 mesh screen to substitute the 60 μm particles in numerical simulation.

The volume of 120 μm particles was eight times that of 60 μm particles. The volume of 120 μm particles put into the hydrocyclone was eight times that of 60 μm particles to ensure an equal number of both particles.

Table 5 shows the standard volume concentration adopted. The concentrations for 60 and 120 μm particles were 5 kg·m−<sup>3</sup> and 40 kg·m<sup>−</sup>3, respectively.


**Table 5.** Experiments for particle separation.

**Experiment 1.** *Sand particles (61–106 μm; 9 kg) were screened out using a 150 to 240 mesh screen to substitute the 60 to 106 μm irregular particles in numerical simulation, and a separation simulation was performed on them.*

**Experiment 2.** *Three experiments were carried out, including two independent single separations of 60 and 120 μm particles and one mixed separation of both particles.*

### **3. Analysis of Experimental Results**

### *3.1. Flow Field, Flow Rate and Particle Distribution*

The separation process of water and sand in a hydrocyclone is as follows: the turbulent flow of water drives the sediment particles to move, as a result of which the latter eventually flow out from the upstream and downstream separately. It can be seen that turbulent flow of water plays an important role in the separation of sediment particles. After it flows in from the inlet, the water flows into a swirling turbulent flow in the cavity, and finally carries the sediment particles out of the upstream or downstream separately.

Figure 4 shows the velocity of water flow in the hydrocyclone along the *X*-, *Y*- and *Z*-axes, the resultant velocity of water flow (scalar quantity obtained by calculating the values of three coordinate axis), and the simulated result of water flow pressure distribution.

**Figure 4.** Structure of the cyclone, water flow velocity (m·s−1) and pressure (Pa) distribution: (**a**) structure of hydrocyclone; (**b**) flow velocity in X-axis; (**c**) flow velocity in Y-axis; (**d**) flow velocity in Z-axis; (**e**) water flow velocity distribution; (**f**) water flow pressure distribution.

Figure 4a shows the structure of the hydrocyclone, Figure 4b, 4c and 4d are the water flow velocity (m/s) distributions along the *X*-, *Y*- and *Z*-axes, respectively. Figure 4e and 4f presents the water flow velocity distribution and the water flow pressure distribution, separately. The velocity distribution along the *X*-axis is symmetrically distributed on both sides centered with the *Z* axis. The velocity distribution along the *Y*-axis exhibits periodic positive and negative phases. With regard to the speed distribution along the *Z*-axis, the velocity near the wall of the cylindrical section is low, and particles there flow from the downstream. The water flows through a zero-speed buffer zone when approaching the axis center, and eventually flows to the upstream. The general trend is that the water flow velocity in the cavity changes significantly at the intersection between cylinders, and the flow velocity increases near the downstream (as shown in Figure 4e). The pressure distribution in the cavity did not show obvious changes, but changes significantly when near the outlet. By observing the water flow velocity and pressure distribution, it can be known that the water flow carries the sediment particles from the inlet, and being blocked by the wall, particles rotate in circles in the cyclone chamber, and finally flow out from the upstream or the downstream. Since the water carrying capacity is influenced by the shape and size of the particles themselves, the hydrocyclone is usually used for sediment particle separation.

Figure 5 shows the position distribution of the particles in the hydrocyclone during the simulation. Figure 5a shows the separation simulation of 60 and 120 μm particles, and Figure 5b shows the simulation of the effect of the cross-section change of sand particles in the separation. Figure 5c shows the distribution of sand particles along the *Z*-axis and Figure 5d shows the particle distribution at the D-D section in Figure 5b.

**Figure 5.** Sand particle distribution in simulation: (**a**) the separation simulation of 60 and 120 μm particles; (**b**) the simulation of the effect of the cross-section change of sand particles in the separation; (**c**) the distribution of sand particles along the Z-axis; (**d**) the particle distribution at the D-D section in (**b**).

From Figure 5c,d, the overall distribution of sand particle in the cylindrical section and part of the conical section can be obtained, and the closer of particles are to the wall, the higher the density is. The radial (*X*-axis) concentration distribution, as well as the velocity and the separation results in Figure 5b are shown in Tables 6 and 7, respectively. The axial velocity of the *Z*-axis, the radial (*X*-axis) distribution percentage, and the separation results in Figure 5a are shown in Figure 6, Tables 8 and 9, respectively.

### *3.2. Relationship between Maximum Projected Area and Separation Results*

### 3.2.1. Radial Concentration Distribution of Particles

Table 6 shows the concentration distribution of different sand particles in the radial area of Figure 1. The general trend of sand particle distribution was that the higher the concentration distribution and the nearer the sands to the wall, the more concentrated the large sand particles, especially the 106 μm sand particles, which accounted for 86.93% of the total area.


**Table 6.** Concentration distribution of different sand particles in the cylindrical section.

The volume of 60 μm particles was the same as that of particle types C, B, and A. Their concentrations were 69.73%, 65.44%, 57.58%, and 51.43%, respectively, in the area with 4 to 5 mm radial distance near the wall. Therefore, the concentration of sand particles with the same volume tended to decrease gradually near the wall of the hydrocyclone with the increase in the maximum projected area.

The concentrations of 60, 70, 90, and 106 μm particles in the area with 4 to 5 mm radial distance were 69.73%, 76.00%, 82.62%, and 86.93%, respectively. That is to say, the concentration of spherical particles tended to increase in the direction of the hydrocyclone wall with the increase in the volume.

Taking the 2-3 area as the center, the concentrations of types C, B, and A in the area with 0 to 2 mm radial distance were 4.01%, 9.79%, and 14.14% higher than those of 70, 90, and 106 μm particles, respectively. Furthermore, their concentrations in the area with 3 to 5 mm radial distance were 5.97%, 15.00%, and 23.80% lower, respectively. Therefore, for the particles with the same volume and maximum projected area, the higher the concentration in the center of the hydrocyclone, the more likely the central part (the upwelling) flowing out from the upstream; for the particles with the same maximum projected area, the bigger the volume, the higher the concentration near the hydrocyclone wall, and the closer to the side wall, the more likely they might flow out from the downstream.

### 3.2.2. Separation Results in Numerical Simulation

Table 7 shows the separation results of the upstream and downstream with changes in the volume and maximum area.


**Table 7.** Separation results of particles from the upstream.

tupstream, the average residence time of particles crossing the upstream; tdownstream is the average residence time of particles crossing the downstream; Vupstream is the average absolute velocity of particles crossing the upstream; Vdownstream is the average absolute velocity of particles crossing the downstream; Lupstream is the average path length of particles crossing the upstream; Ldownstream is the average path length of particles crossing the downstream; and *p* is the percentage of the particles crossing the upstream in the separated particles, which was obtained using Equation (19).

As shown in Table 7, 9.4%, 7.4%, 3.2%, and 0.9% of spherical particles with diameters of 60, 70, 90, and 106 μm flowed out from the upstream, respectively. Therefore, it was harder for larger spherical particles to be separated from the upstream. The 60 μm particles had the same volume as types C, B, and A particles. As the maximum projected area increased gradually, 9.4%, 14.9%, 17.8%, and 24% of them flowed out from the upstream, accordingly. Thus, it was easier for the particles with the same volume to be separated from the upstream as the maximum projected area increased.

The 60 μm particles and types C, B, and A particles flowed out from the upstream in less time and traveled a shorter distance with the increase in the maximum projected area. For instance, the particles of type A took 0.09 s less compared with 60 μm particles and the average distance decrreased by 0.03 m.

However, the 60, 70, 90, and 106 μm particles needed more time to flow out from the upstream and travel a longer distance with the increase in their size. Among them, the 106 μm particles took the most time, 0.09 s more compared with 60 μm particles on average, and the average distance increased by 0.03 m.

For particles obtained from the upstream, the average velocity of 60 μm spherical particles and types A, B, and C particles was greater than that of 70, 90, and 106 μm particles. However, the average velocity of particles obtained from the downstream was in contrast to those from the upstream. The absolute velocity of different particles was less than that of water flow. Therefore, the smaller the particles or the larger the projected area, the stronger the ability to flow.

### *3.3. Separation Mechanism of Single and Mixed Particles*

### 3.3.1. Distribution of Axial Velocity and Radial Concentration

Figure 6 reflects the radial velocity of 60 and 120 μm particles. Three curves are seen from the top downward. The first one is about the average velocity of 120 μm particles in a positive direction, the direction of downstream. (The average velocity of 120 μm particles in single and mixed separations is shown by the same curve as it is the same in different separations, and 60 μm particles have no influence on 120 μm particles.) The second curve shows the average velocity of 60 μm particles in mixed separation in a positive direction (the direction of downstream). The third one is about the average velocity of 60 μm particles in single separation in a positive direction of *Z* axis. Figure 6b shows the partially enlarged image of Figure 6a (the former is five times bigger than the latter).

**Figure 6.** Average velocity distribution in a positive direction of *Z* axis: (**a**) radial velocity of 60 and 120 μm particles; (**b**) partially enlarged image of (**a**).

The velocity and direction of water flow significantly changed at 25–81 mm in Figure 6, which is the transition from a cylinder to a cone. Particles slowed down in a positive direction on the *Z* axis due to a lower water velocity, and the average velocity of different-sized particles also decreased. Thus, particles in this section mostly flow slowly in the whole separation process.

The velocity of 60 and 120 μm particles on *Z* axis had four changing stages: rising, falling, flattening, and rising again. Compared with a single separation, the average velocity of 60 μm particles increased by 25.74% at 25–81 mm in mixed separation, which was quite obvious at 25–51 mm. However, the velocity tended to be consistent at 0–25 mm and 81–151 mm.

The velocity was the lowest at 25–81 mm in the whole process. As shown in Figure 1, it was the border of cylinder and cone sections, 25–51 mm in the cylinder section and the rest in the cone section.

Table 8 shows the distribution of particle concentration at A to B in Figure 4. It was found that 93.23% of 120 μm particles in single separation were concentrated 4–5 mm away from the center of the circle in the radial direction. The percentage of particles nearest to the center of the circle was 1.33%. However, 69.73% of 60 μm particles in single separation were concentrated at 4–5 mm, and the percentage of particles nearest to the center of the circles was 2.65%. When the particles of these two sizes were mixed, the radial concentration distribution of 120 μm particles was in accordance with that of single separation generally. However, the radial concentration distribution of 60 μm particles in single separation was quite different from that of mixed separation. The difference was reflected in the section of 4–5 mm, and the percentage reached 20.73%.

In other words, influenced by 120 μm particles, 20.73% of 60 μm particles were concentrated in the aforementioned section in mixed separation.


**Table 8.** Concentration distribution of 60 and 120 μm particles in the cylinder section.

### 3.3.2. Particle Separation Results

The calculation of particle separation took 4 s (see Table 9). The number of particles that flowed into the hydrocyclone was more than the sum of particles separated from upstream and desilting port in 0–3 s; the separation was unstable. However, it reached a dynamic balance because the number of particles flowing into the hydrocyclone was close to that of the particles separated from it in 3–4 s. When 60 μm particles were mixed with others, the average time spent in separation from the downstream was 0.2 s less than that in single separation while the separation increased by 4.1%. Furthermore, 120 μm particles were not separated in either single separation or mixed separation from the upstream.

**Table 9.** Separation of 60 and 120 μm particles.


X and D, the number of X (60 μm particles) and D (120 μm particles) was 1000 per second; tavearge is the average time duration in which particles were separated from the downstream staying in the hydrocyclone in 3–4 s; Ndownstream is the number of particles passing the downstream; Nupstream is the number of particles passing the upstream in 3–4 s; and *p* is the percentage of particles in the downstream, which was calculated using Formula (19).

### **4. Experimental Tests**

### *4.1. Test on the Effect of Cross-Section Change on Separation*

Figure 7 depicts the observation results of sand samples in the experiment. Figure 7a shows the representative sand particles screened using a vibrating screen. Both particles 1 and 2 could pass the 110 μm hydrocyclone screen, but they were not typical spherical particles, and their volume was less than that of spherical particles of the same size. Figure 7b shows the typical comparison of thickness. Obviously, particles 1 and 3 were thicker than particle 2 in terms of light transmittance. Figure 7c shows the particles with different projected area and thickness. Among them, the volumes of particles 1, 2, 3, and 4 might be the same. Figure 7d shows statistics on the maximum projected area of typical lamelliform particles separated from the upstream. The maximum projected area of particles 1, 2, 3, and 4 in Figure 4 was 7968.39, 8346.47, 8596.55, 3826.37, and 3949.87 μm2, respectively. Accordingly, the size of spherical particles was 100.75, 103.11, 104.65, 69.81, and 70.93 μm, whose volume was less than that of spherical particles with the same projected area. Smaller particles were more than larger particles. Moreover, most were lamelliform particles with a similar size.

**Figure 7.** Distribution of particle characteristics (×1 k).

Different particles (10 groups) were observed with a transflective polarized microscope whose number was 59XC-PC before the experiment. It was found that most of 500 particles had a bad light transmittance and only a small number of them (10–30 particles) had a good light transmittance. However, about 160–210 among 500 particles obtained from the upstream were good in light transmittance (because they were relatively thin) after the experiment.

The results indicated that the larger the size of particles, the less the particles were separated from the upstream. As shown in Figure 8, the percentage of 60 and 106 μm particles was 27.75 and 1.91%, respectively.

Therefore, smaller and thinner particles were separated from the upstream more easily. In other words, the process of separating sand particles in the hydrocyclone led to the classification of particles of different sizes, which was equivalent to the filtering of particles of different shapes (lamelliform particles). The experimental result was consistent with the regularities reflected in numerical simulation.

**Figure 8.** Percentage of particles with different sizes on the upstream.

### *4.2. Separation Test of Sand Particles 60 and 120 μm*

The concentration of 60 and 120 <sup>μ</sup>m sand particles was set as 5 kg·m−<sup>3</sup> and 40 kg·m<sup>−</sup>3, respectively, to ensure an equal number of sand particles put into the hydrocyclone. The capacity of the container was 1 m3. The separation experiment lasted for 16 min.

Table 10 shows the results: 120 μm sand particles were separated from the downstream in both single and mixed separations. However, 312 g of 60 μm sand particles were separated from the upstream in single separation, and 255 g were separated in mixed separation. It was concluded that larger particles improved the separation of smaller ones from the upstream in mixed separation, which was consistent with the regularities reflected in numerical simulation.


**Table 10.** Separation experiments for 60 and 120 μm sand particles.

### **5. Discussion**

### *5.1. Influence of Maximum Projected Area and Volume Changes on Separation Results*

Several parameters affect the performance of hydrocyclones, including the inlet velocity, solid content, liquid-phase viscosity, particle size and shape, and so forth. Saber et al. [20] believed that the separation effect would improve with the increase in the size of spherical particles, which was obviously better than that of plate-like particles. Also, the separation effect declined with the increase in the size of plate-like particles with diameters more than 60 μm. Abdollahzadeh et al. [21] thought that the separation effect decreased along with the increase in flatness when the particle diameter was more than 10 μm. Endoh et al. [22] found that particles with a large size and high proportion were separated from the upstream. Wills et al. [23] held the view that plate-like particles, such as mica, were often discharged upstream even though they were relatively coarse. Kashiwaya et al. [17] believed that the shape of particles would influence the drag coefficient and separation effect. The drag coefficient would increase when the particles became flatter. Then, the increased drag force would increase the influence of water movement on particles. Wang et al. [24] found that the change in particle velocity fluctuated evidently as the size of particles became smaller and the drag force increased. Some particles passed through the locus of zero vertical velocity, joined the upward flow, and then moved out from the upstream. In this study, particles of types C, B, and A, corresponding to 70, 90, and 106 μm particles, respectively, could enter the upward flow and then flow out of the upstream easily as particles became smaller and drag force increased. Therefore, the corresponding percentage dropped sharply in the area with 3 to 5 mm radial distance and increased in the area with 1 to 3 mm radial distance.

Drag force is proportional to the maximum projected area of particles along the direction of fluid flow (Formula (19)), which means that the drag force is strengthened with the increase in the maximum projected area of particles. Hence, the maximum projected area of particles has a certain influence on the separation results of the hydrocyclone. Kashiwaya et al. considered that when Reynolds number reached a certain level (above 50), the maximum projected area of particles was a major influencing factor. They also claimed that the radial variation in the tangential velocity in the hydrocyclone gave rise to a moment on plate-like particles and the particles settled, maintaining the orientation as their broad planes became perpendicular to the radial direction. In this study, the pressure gradient force in the upwelling direction toward the center of the hydrocyclone in the same position increased with the increase in the cross-sectional area. Also, the particles could be easily pushed into the upwelling and then out of the upstream. Figure 9 shows that the volume of 60 μm spherical particles was the same as that of particles of types C, B, and A, but the cross-sectional area was larger and hence the track was simpler in the hydrocyclone.

Figure 9 shows the representative motion tracks of sand particles with different cross-sectional areas when they moved from the upstream. The figure shows that 60 μm spherical particles had the same volume as that of particles of types C, B, and A. With the increase in the maximum projected area, the tracks became simpler and the time needed for separation became shorter. Furthermore, 60, 70, 90, and 106 μm particles were all spherical. The tracks separated from the upstream become more complex with the increase in the volume of particles. That is, it was harder for particles to separate from the upstream as the time taken for separation increased. Tracks (2) and (8), (3) and (7), and (4) and (6) belong to the particles with the same maximum projected area. The smaller the volume of particles, the simpler the tracks when separated from the upstream, indicating that the time taken for separation was shorter, and it was harder for particles to separate from the upstream. Hence, the maximum projected area might have an obvious and a stable influence on the motion tracks of particles, hence affecting the separation effect.

**Figure 9.** Representative motion tracks of sand particles.

### *5.2. Analysis of the Following Phenomenon of Fine Particles*

Zhang et al. [19] believed that the movement of particles was controlled by the combined effect of the radial fluid drag force, pressure gradient force, and centrifugal force. Of these, the fluid drag force of the inward direction increased exponentially with decreasing size of small particles. A part of the small particles was pushed into the center of the hydrocyclone and then flowed out of upstream. Xu et al. [25] held the view that large particles moved toward the wall under the effect of inertial force and gravity. In this test, a part of 60 μm particles flowed out of the upstream, while the 120 μm particles did not. However, the following phenomenon enhanced the density of the 60 μm particles near the wall, where most 120 μm sand particles were present, thereby reducing the possibility of 60 μm particles flowing into the center of the hydrocyclone and out of the upstream. Therefore, the force on particles in the hydrocyclone was complicated. Zhu et al. [26] noted that the comprehensive effect of the fluid drag force, pressure gradient force, virtual mass force, and gravity on particles should be taken into consideration. Peng et al. used the lattice Boltzmann method to successfully simulate the shallow water flows and river bed evolution [27–34] as well as two-phase flows [35,36].

Tsuji et al. [37] believed that the large leading particles had a huge drag force on small trailing particles, which had little or even no drag force on leading particles in turn. The experiment showed that 120 μm particles had a great impact on 60 μm particles in terms of distribution, velocity, or separation effect. However, no 120 μm particles separated from the upstream. Furthermore, 120 μm sand particles could influence 60 μm particles to a great degree, but the latter's influence on the former could be neglected, which was consistent with the findings of other studies. Meanwhile, Zhu et al. [38] claimed that Re could affect drag, and the increase in Re could obviously weaken the influence of leading particles on small ones. Zhang et al. [18] held the view that the following drag of small particles toward large ones would decrease with the increasing distance, and it would be affected by particle's ratio and the local wake flow compared with that between the particles and the incoming fluid stream. Schubert et al. [39] considered that this phenomenon was due to particle interaction, while the finer particles were captured by the larger particles and the effect was more prominent with the increase in the fraction of large particles. In this study, Re did not change at all under the force of inlet pressure, which was the only one used in this experiment. However, the following phenomenon was quite obvious from the perspectives of process and results.

Under the ideal circumstances, the deposition velocity of spherical particles is proportional to particle diameter (Formula (13)). The velocity is high with the increase in the particle diameter. The numerical calculation and experimental results in this study showed that larger particles were faster than smaller particles in terms of deposition velocity in the hydrocyclone. CFD–DEM coupling was adopted to capture the following phenomena of the movement of small particles toward large ones (see Figure 10). Figure 10a,b presents the motion tracks of 120 and 60 μm particles in the downstream in mixed separation, whereas Figure 10c is the enlarged drawing of the complex segment in Figure 10b. In the figure, among particles around the transition section from a cylinder to a cone, 60 μm particles firstly entered the hydrocyclone, and their motion track was chaotic. After meeting with the 120 μm particles, the 60 μm particles changed their direction and were separated from the downstream together with the 120 μm ones. This phenomenon also explains why the average velocity of mixed separation of 60 and 120 μm particles near the boundary of cylinder and cone was higher than that of single separation (Figure 6).

**Figure 10.** Representative motion tracks of particles.

### **6. Conclusions**

In this study, numerical simulation and experimental methods were applied to reveal the separation mechanism of the hydrocyclone for different-sized particles and the influence of maximum projected area on the separation effect. The conclusions were as follows:

The percentage of 60, 70, 90, 106, and 120 μm spherical particles separated from the upstream was 9.4%, 7.4%, 3.2%, 0.9%, and 0%, respectively, indicating that spherical particles were more easily separated from the downstream as their size increased.

For 60 μm spherical particles and types C, B, and A particles with the same volume, the separation percentages from the upstream were 9.4%, 14.9%, 17.8%, and 24.0%, respectively, implying that the running space was closer to the upward flow in the hydrocyclone as the maximum projected area increased. Hence, particles with the same volume were easier to be separated from the upstream.

The axial velocity of 60 μm particles mixed with 120 μm particles in separation increased by 25.74% compared with that in the single separation of 60 μm particles near the transition section from

a cylinder to a cone. The concentration near the wall, where mainly 120 μm particles were present, increased by 20.73%, and the separation effect in the downstream improved by 4.1%. Therefore, large particles could optimize the separation effect of small particles from the downstream.

**Author Contributions:** Software, Z.T.; methodology, N.L. and F.W.; writing—original draft preparation, Z.T.; writing—review and editing, L.Y.; funding acquisition, L.Y., L.C. and N.C.

**Funding:** This research was funded by the National Natural Science Foundation of China (NSFC) (No. 51769009, 51379024,51809022) and Funding Program for University Engineering Research Center of the Yunnan Province in China.

**Acknowledgments:** The authors appreciate the financial support received from the National Natural Science Foundation of China (NSFC) (No. 51769009, 51379024, 51809022) and Funding Program for University Engineering Research Center of the Yunnan Province in China.

**Conflicts of Interest:** Declare conflicts of interest or state.

### **Nomenclature**


### Greek letters


### **References**


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

### *Case Report* **A Comprehensive Method of Calculating Maximum Bridge Scour Depth**

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


Received: 4 October 2018; Accepted: 1 November 2018; Published: 3 November 2018

**Abstract:** Recently, the issues of scour around a bridge have become prominent because of the recurrent occurrence of extreme weather events. Thus, a bridge must be designed with the appropriate protection measures to prevent failure due to scour for the high flows to which it may be subjected during such extreme weather events. However, the current scour depth estimation by several recommended equations shows inaccurate results in high flow. One possible reason is that the current scour equations are based on experiments using free-surface flow even though extreme flood events can cause bridge overtopping flow in combination with submerged orifice flow. Another possible reason is that the current practice for the maximum scour depth ignores the interaction between different types of scour, local and contraction scour, when in fact these processes occur simultaneously. In this paper, laboratory experiments were carried out in a compound shape channel using a scaled down bridge model under different flow conditions (free, submerged orifice, and overtopping flow). Based on the findings from laboratory experiments coupled with widely used empirical scour estimation methods, a comprehensive way of predicting maximum scour depth is suggested which overcomes the problem regarding separate estimation of different scour depths and the interaction of different scour components. Furthermore, the effect of the existence of a pier bent (located close to the abutment) on the maximum scour depth was also investigated during the analysis. The results show that the location of maximum scour depth is independent of the presence of the pier bent but the amount of the maximum scour depth is relatively higher due to the discharge redistribution when the pier bent is absent rather than present.

**Keywords:** bridge scour; sediment transport; submerged flow; physical hydraulic modeling

### **1. Introduction**

When a bridge is constructed in a river, the flow pattern around the bridge changes because a unique flow field develops locally around the bridge pier and abutment. Furthermore, reduced flow area through the bridge opening by the existence of embankment/abutments at both and/or one side of the river results in higher flow velocity due to the acceleration. This unique flow field with the higher velocity can seriously damage bridge foundations. Thus, if the depth of the foundation is not deep enough, the chance of bridge failure becomes higher.

A bridge can fail because of several causes, such as earthquake, wind, and flooding. Among them, bridge scour is the biggest reason of bridge failure [1,2]. For example, about 60% of bridge failures of total bridge collapse in the United States since 1950 have been related to the scour of bridge foundations [3]. The Colorado Department of Transportation (CDOT) estimated a minimum of 30 state highway bridges were destroyed and twenty were seriously damaged by floods in the year of 2013 [4]. In Nepal, due to the degradation of bed materials during the 2014 flooding, the foundation of the

highway bridge over the Tinau River was seriously exposed [5]. As explained in the above examples, it is justified to say that bridge scour is one of the main bridge safety problems all over the world. Thus, accurate prediction of scour at the bridge foundation becomes the primary aim of engineers for the safety of a bridge.

Since the late 1950s, numerous studies on scour around bridges have been conducted and have generated formulas for equilibrium scour depth estimation at bridge foundations [6–9]. Although using equilibrium scour depth around a bridge foundation is a reasonable practice to design a bridge, sometimes it shows an overly conservative design compared to the field measurements. Contrary to the conservative design, a study in South Carolina found that the observed scour depth was greater than the equilibrium scour depth based on using 100-year flooding even if there had not been a 100-year flooding since the bridge was built [10,11]. Repeated occurrences of smaller flooding events might cause scour that was greater than the scour estimation in South Carolina. Moreover, some researchers have used computational fluid dynamics (CFD) simulations and found that the scour depth was under-predicted by the steady-state calculations while it was over-predicted by the unsteady-state [12,13].

One of the possible reasons of lacking an accurate scour prediction method is that many of the scour predictor equations were derived from simplified laboratory studies under free flow cases. The extreme amount of water associated with huge flood events can result in a complex flow field around the bridge under bridge overtopping flow in combination with submerged orifice flow. In addition to the complex flow field, the irregular shape of river geometry as well as non-uniform sediment distribution cannot be reproduced in a simplified-idealized laboratory setting such as in rectangular flume as in the previous research. Also, most of the scour predicting equations are 2nd or 3rd order which may not accurately predict scour depth as the scouring process is a complicated scholastic phenomenon.

Another possible reason for the inaccurate prediction of scour depth is that the current practice of total scour depth assumes that contraction and local scour are independent processes. Contraction scour is the consequence of flow acceleration due to contraction in the flow area, while local scour is caused by local vortex structures around the base of the obstruction. However, when the bridges are constructed in the river, these two-flow patterns tend to occur concurrently, which make local scour and contraction scour time dependent [14]. Thus, to predict maximum scour depth around a bridge foundation, a single equation should be developed rather than two separate equations for different types of scour.

Hence, the main objective of this study is to develop a single equation that can be used to predict maximum scour depth where different types of scour occur simultaneously. A scale down laboratory experiment was conducted to find the interaction between different scour components and the result was applied to and compared with the most widely used scour equations in the US (Colorado State University (CSU) and Melville-Sheppard (M/S) equations) to suggest an improved way of calculating the maximum scour depth. These equations predict maximum pier scour depths. Basic applications include simple pier substructure configurations and riverine flow situations in alluvial sand-bed channels. The CSU equation has been used for bridge scour evaluations and bridge design for countless bridges in the U.S. and worldwide. The M/S equation combines pier geometry, shape, and angle of attack to compute an effective pier width, *a*\*, and also distinguishes between clear-water and live-bed flow conditions.

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

### *2.1. Physical River Modeling*

As shown in Figure 1, a 1:60 scale down physical model was constructed in the hydraulics laboratory including full river bathymetry of the Towaliga River bridge at Macon, Georgia by using Froude number similarities between the model and the prototype [15–17]. A relatively small bridge was selected in this experiment because a large number of smaller bridges can fail during extreme hydrologic events. The drainage area of the selected site at the bridge is 816 km2. Discharge was estimated as 1700 m3/s by the U.S. Geological Survey for Tropical Storm Alberto in 1994 which is larger than 500 year flooding. During this historic event, severe overtopping of the bridge and significant scour around the pier bents in the left floodplain occurred, but the embankment and bridge remained sufficiently intact until repairs could be made to the scoured area. The actual scour depth measured during the storm events can be found in Hong and Sturm [16,17]. As shown in Figure 1, each of the pier bents consists of two in-line rectangular columns having a width of 1 m, and they were modeled inclusive of pile caps and piles. A bridge deck width of 13 m, in accordance with standard two-lane roads was also constructed in the laboratory.

**Figure 1.** Towaliga River bridge in the field and model in the laboratory.

The approach section of the bridge was 7.3 m long followed by a 6.09 m long moveable working bed section. The approach section was filled with uniform size of small gravel (*d*<sup>50</sup> = 3.3 mm) to make a fully rough turbulent approach flow. In the moveable working bed section, the full depth was filled with sand with median diameter size (*d*50) of 0.53 mm. The bridge model was constructed within the moveable working bed section for the scour experiment.

For selecting sediment size in the laboratory study, recently developed scour modeling methodology was applied [2,18–20]. Selecting sediment size in the laboratory was mostly related to the ratio of pier size to sediment size (*a*/*d*50). In the field, the effect of sediment size, *d*50, has not been considered to be important to pier scour because of very large values of *a*/*d*<sup>50</sup> [21]. However, in the laboratory, pier scour depth to pier size ratio (*dps*/*a*) tends to increase with *a*/*d*<sup>50</sup> up to a maximum at *a*/*d*<sup>50</sup> ≈ 25 and seemingly becomes independent of the ratio when *a*/*d*<sup>50</sup> is greater than 50 [19]. However, another researcher has suggested that *dps*/*a* may decrease at very large values of *a*/*d*<sup>50</sup> based on experiments in a large flume [22]. Therefore, the sediment size, *d*50, was chosen such that the ratio of pier size to sediment size, *a*/*d*<sup>50</sup> was in the range of 25–50 where it has negligible influence on pier scour. Furthermore, the ratio of the approach velocity to the critical velocity which concludes the condition of clear water scour regime was also an important factor to choose sediment size in the experiment. Considering all the above conditions, the sediment size for this experiment was chosen as *d*<sup>50</sup> = 0.53 mm with *σ<sup>g</sup>* = 1.2, where *σ<sup>g</sup>* is geometric standard deviation of the particle size distribution.

Within the working moveable bed section, the erodible embankment/abutment was constructed at both sides of the channels by using erodible fill with rock-riprap protection in order to reproduce the same conditions as in the field [16,23]. Furthermore, the existing roadway and bridge deck were modeled based on the scale-ratio and put on top of the embankments. With this approach, submerged orifice flow with and without overtopping during the extreme amount of flooding as observed in the field could be simulated in the laboratory. As shown in Figure 2, eight pier bents were also constructed through the bridge to support the bridge deck and for the local scour experiment. Each of the pier bents consisted of two in-line rectangular columns. Abutment structures were constructed and buried at both sides of the erodible embankments.

**Figure 2.** Geometry of compound channel for (**a**) plan view with velocity measurement locations; (**b**) cross section view at bridge.

### *2.2. Experimental Procedure*

After completion of the model structure, the flume was slowly filled with water from a downstream supply hose to saturate the sand without disturbing the initial bottom contours. Then bottom elevations were measured in detail throughout the entire working section using an Acoustic Doppler Velocimeter (ADV). After that, a larger flow depth than the required value was set by the tailgate, then discharge was increased slowly to prevent initial scour while setting up the test discharge. Then the tailgate was lowered to achieve the desired depth of flow. In the meantime, a point gauge was used to measure the flow depth to measure the targeted water depth. Once the desired flow rate and water depth were achieved, scouring was continued for 5 to 6 days until equilibrium (change in scour depth less than 2% within 24 h) was reached. After reaching equilibrium condition, the entire bed elevation (bed elev.) was measured by a point gauge and the ADV in detail to obtain accurate contours after scour. After finishing the moveable bed experiment, the surface of the movable bed was fixed by spraying polyurethane. In the fixed bed conditions, the velocities were measured by ADV in the approach section and bridge upstream (U/S) and downstream (D/S) section. During the velocity measurements, typical correlation values in the experiments were greater than 90% and the Signal Noise Ratio (SNR) was greater than 15. The sampling frequency of the ADV was chosen to be 25 Hz with a sampling duration of 2 min at each measuring location. More detailed measurements techniques using ADV can be found in several other articles [24–29].

### *2.3. Assessment of Reference Scour Depth*

Two well established theoretical pier scour equations were used to decide the reference scour depth at the pier for this research. One of the most commonly used pier scour equation in the United State is the CSU equation (also known as the HEC-18 equation). The CSU equation includes a correction factor for pier shape, angle of attack of flow, and the bed conditions. The CSU equation was initially developed from a laboratory data set measured by several researchers [6,30,31]. After the initial development, the CSU equation has been progressively modified over the years and is currently recommended by the Federal Highway Administration (FHWA) for estimating equilibrium scour depths at simple piers as follow [6]:

$$\frac{d\_{\rm csu}}{\chi\_2} = 2K\_1 K\_2 K\_3 \left(\frac{a}{\chi\_2}\right)^{0.65} Fr\_2^{0.43} \tag{1}$$

where, *K*1, *K*2, *K*<sup>3</sup> are the correction factors for pier nose shapers, for angle of attack of flow, and for bed conditions, respectively; *Fr*<sup>2</sup> = Froude number at the pier face = *V*2/(*gY*2) 0.5; *a* = Pier width; *dcsu* = Equilibrium pier scour depth; *Y*<sup>2</sup> = Flow depth directly upstream of the pier; *V*<sup>2</sup> = Mean velocity of flow directly upstream of the pier; *g* = Acceleration of gravity. Several evaluations of the CSU equation through laboratory and field datasets have shown this equation may perform better than other bridge predictive equations [32–34]

Another commonly used pier scour prediction method is the Melville–Sheppard or M/S equation [6,34,35]. A new variable is introduced in the M/S equation, effective pier width, *a*\*, which shows the combined effect of pier geometry, shape, and angle of attack. The equation is as follows:

$$\frac{d\_{\rm ms}}{a^\*} = 2.5 f\_1 f\_2 f\_3 \tag{2}$$

where, *f* 1, *f* 2, *f* <sup>3</sup> are the factors for flow-structure interactions, for flow-sediment interactions, and for sediment-structure interactions, respectively. Here,

$$f\_1 = \tanh\left[\left(\frac{\chi\_2}{a^\*}\right)^{0.4}\right] \tag{3}$$

$$f\_2 = \left[1 - 1.2\left[\ln\left(\frac{V\_2}{V\_c}\right)\right]^2\right] \tag{4}$$

$$f\_3 = \left[\frac{\left(\frac{\mu^\*}{d\_{50}}\right)^{1.13}}{10.6 + 0.4\left(\frac{\mu^\*}{d\_{50}}\right)^{1.33}}\right] \tag{5}$$

where, *V*<sup>c</sup> is the critical velocity for movement of *d*50, calculated by Keulegan's equation.

Using the provided theoretical pier scour equation, pier scour depth was predicted with the variables measured in the experiments. Because theoretical equations were derived from a simple laboratory set-up in a rectangular flume, the result can be used as a reference pier scour depth without any effect of flow contraction. Then, the reference pier scour depth was compared with the measured maximum scour depth around the pier where the flow contraction and local vortex structure interact and lead to the maximum scour depth. The results showed the effect of flow contraction on the pier scour and were used to suggest precise estimation of pier scour depth [8].

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

A total of eight experiments were conducted for this research and the experimental conditions are summarized in Table 1 where *Q* is the total discharge and *q*2/*q*<sup>1</sup> is the discharge contraction ratio which can be used as a key-independent variable that accounts for flow redistribution and resulting flow acceleration through a bridge section. *V*1, *Y*1, and *V*2, *Y*<sup>2</sup> are cross-sectional mean velocity, and depth of floodplain in the approach and at the upstream face of the bridge section, respectively. Also, maximum scour depth expressed as flow depth measured at the location of maximum scour (*Y*m) and longitudinal (flow direction) distance measured from the upstream face of the bridge to the location of maximum scour depth (*X*) (see in Figure 2a) are presented in Table 1. A better idea of the scour distribution at the bridge cross section in these experiments can be seen in Figure 6 as a reference. Run 9 and 10 were selected from the previous study [25] and used to validate the results in this research so that the suggested method can be used for different hydraulic and geometry conditions. As shown in Table 1, runs 1, 4, 5, and 6 were performed in free (F) flow condition; whereas, the remaining runs were conducted in the submerged flow conditions with overtopping (OT) and without overtopping (SO).

Close to the toe of the abutment where the abutment side slope ends and the first pier is introduced, these locations are the most vulnerable to scour because abutment scour, pier scour, and contraction scour occur simultaneously. Thus, to find the effect of pier on the maximum scour depth, pier bent #7

was removed from the river model and experiments were conducted in runs 7 and 8 with the exact same flow conditions as in runs 2 and 3. Even if the velocities and scour depths were measured in the entire working moveable bed section, only the floodplain flow variables and scour depths are presented in this paper because the maximum scour depth occurred on the floodplain in all our cases.

For the velocity measurements at each cross-section, point velocities were measured along multiple vertical transects through the entire cross-section, and at each vertical transect. Point velocities were measured at minimum of four points and maximum ten points vertically depending on the depth of water. Then, the depth-averaged velocity was determined by the best fit of the logarithmic velocity profile in each vertical transect. Based on the depth-averaged velocity, the cross-sectional mean velocity of the floodplain was calculated for the approach and bridge section as presented in Table 1.


**Table 1.** Experimental conditions and location of maximum scour depth measurement.

Notes: F = free flow, SO = Submerged orifice flow, OT = Overtopping.

### *3.1. Measurement of the Maximum Scour Depths*

Initial bottom elevations were measured throughout the test (moveable bed) section before the experiment and then the final bottom elevations were measured at the same locations in equilibrium conditions. Then, to find the location and magnitude of the maximum scour depth, bed elevations after and before the scouring were compared. The maximum scour depth in all of the laboratory experimental runs was found within the bridge and close to the bridge pier (Table 1).

### *3.2. Prediction of Maximum Scour Depth*

The main assumption in this research is that maximum scour depth consists of theoretical pier scour depth and additional scour due to flow contraction. Based on this assumption, the maximum scour depth can be calculated by the equation,

$$\text{Max. Score depth} = \text{Theoretical pier score} + \text{Additional score by flow contraction} \tag{6}$$

The schematic diagram is given in Figure 3 to define the variables used for the analysis. The theoretical pier scour depth can be decided using the CSU (*dcsu*) or M/S equation *(dms*) with measured flow variables and the results are shown in Table 2. Then, the "Additional scour by flow contraction" can be determined by measured the water depth at the deepest point and theoretical pier scour depth as follows;

$$Y\_{m \cdot csu} = Y\_m - d\_{csu} \tag{7}$$

$$Y\_{m\text{-ms}} = Y\_m - d\_{ms} \tag{8}$$

where, *Ym-csu* and *Ym-ms* are the calculated water depths at the deepest location subtracted from the total water depth at the point by the theoretical pier scour depth of CSU (*dcsu*) or M/S equation (*dms*) respectively, which stands for the resulting water depth due to the additional scour. The variables

required for calculation of the reference pier scour depth using the CSU and M/S equations are listed in Tables 1 and 2. Here, *θ* is the flow angle of attack at the pier face which is required to estimate correction factor, *K*<sup>2</sup> and effective pier width, *a*\* for the CSU and M/S equations, respectively. The additional scour components are also shown in Table 2 in terms of non-dimensional variables.

**Figure 3.** Schematic diagram for calculation of maximum scour depth.


**Table 2.** Summary of experimental results to calculate maximum scour depth.

If the assumption is correct, the effect of additional scour due to flow contraction can be represented by the values of *Ym-csu* and *Ym-ms*. Thus, the non-dimensional value of additional scour components calculated using the CSU and M/S equation are compared with the measured discharge contraction ratio in Figure 4a,b respectively. As the value of flow contraction ratio increases, the normalized value of additional scour depth gradually increases. The results clearly reveal that the effect of flow contraction on the additional scour term becomes higher as the value of *q*2/*q*<sup>1</sup> increases. Runs 9 and 10 were also plotted in Figure 4 to validate the findings in different hydraulic and geometry conditions and the result favors the validation as it shows a similar trend of runs 9 and 10 with pressure and free flow data respectively. For both cases using the CSU equation and M/S equation, the pressure flow line has a steeper slope than for the free flow cases and lies above. Because of the vertical flow contraction in addition to the existing lateral flow contraction in pressure flow, the resulting additional scour depth due to the flow contraction is higher than in the free flow cases. Currently, additional experiments are being conducted. With more experimental conditions as well as the data set in Figure 4, the best fit equation will be provided for calculation of scour due to flow contraction under free and pressure flow by a least square regression analysis

It is interesting to note that, the regression analysis results using only with the data set in Figure 4 show higher values of exponents of *q*2/*q*<sup>1</sup> in the case with the M/S pier scour equation for both of the pressure and free flow cases. This finding illustrates that the M/S equation shows a smaller value of pier scour depth in the same flow conditions compared to the CSU equation. The ratios between the theoretical pier scour depths calculated using CSU (*dcsu*) and M/S (*dms*) are plotted with the value of flow intensity factor (*V*2/*V*c) in Figure 5. As shown in Figure 5, the M/S equation shows almost 50% to 75% underestimated values compared to the CSU results. The CSU pier scour equation was developed for live-bed scour and thus ignores the effect of different values of flow intensity and assumes the value is equal to 1; but in the M/S equation, *V*2/*V*<sup>c</sup> is considered as an individual factor to account for the effect of clear water conditions and varies from 0.3 to 1. So, it is obvious that the CSU pier scour equation shows a larger value when applied to the clear water condition. However, there is another important difference between the CSU and M/S equations, which is the consideration of bed material size. However, in our case, the effect of sediment size cannot be considered because all of the experiments are conducted with same sediment size.

**Figure 4.** Effect of flow contraction on additional scour components using (**a**) Colorado State University (CSU) and (**b**) Melville-Sheppard (M/S) equations.

**Figure 5.** Comparison of CSU and M/S pier scour depth in terms of flow intensity.

Finally, a comprehensive procedure for predicting the maximum scour depth with respect to bridge design is introduced:


5. Adding the results from steps 3 and 4 to predict the maximum scour depth for bridge design.

In addition to suggest a comprehensive procedure, the effect of pier bent (located near to the abutment) on maximum scour depth was investigated qualitatively. As shown in Table 1, runs 7 and 8 were conducted with the exact same flow condition as in runs 2 and 3 respectively, but removing pier bent #7 in the bridge section. Figure 6 shows the cross-section comparisons between runs 3 and 8 after scouring. The maximum scour depth occurred at pier bent #6 for both runs 3 and 8. So, the location was independent of the presence of the closest pier bent (#7). However, the amount of maximum scour depth was slightly higher in the case of run 8 where pier bent #7 was absent. A similar trend was observed for runs 2 and 8. Discharge redistribution over the time and interaction between each pier can be the result for this difference of maximum scour depth and needs further study.

**Figure 6.** Comparison of cross-sections for runs 3 and 8.

### **4. Conclusions**

Many investigations have been made attempting to estimate the maximum scour depth and to understand the mechanism of scour around bridge piers. Most of the previous investigations were based on lab experiments using a rectangular channel under free flow. However, due to the recent extreme rainfall events, submerged orifice flow and overtopping flow occur at the bridge frequently, where the flow field around the bridge substructure is more complex than in free flow because of vertical flow contraction in addition to existing lateral flow contraction. Furthermore, most of the natural channel shapes are not rectangular. Also, the current guidelines recommended by HEC-18 assumed that contraction and local scour processes are independent and so they can be determined separately and summed to estimate total scour depth. However, during large flooding events, local scour and contraction scour occur simultaneously and a separate calculation of local scour and contraction scour results in inaccurate scour depth. To overcome the weak points that the current methodology has, laboratory experiments were carried out in a scaled down physical model and a single equation was developed to predict maximum scour depth which can be used without separate calculation of different types of scour components in pressure flow as well as in free flow cases.

Based on the basic assumption of maximum scour depth as summation of theoretical pier scour depth and additional scour depth due to flow contraction, a comprehensive way of predicting maximum scour depth in clear water conditions was suggested. Furthermore, the results demonstrate that the contraction effect on maximum scour depth increases as the flow contraction ratio increases. Also, due to additional vertical contraction in the pressure flow case, the effect of flow contraction on the maximum scour depth shows up larger than in free flow. Another outcome from our investigation concludes that the location of the maximum scour depth is independent of the existence of the closest pier bent but the absence of the closest pier bent increases the scour depth.

Even though this study suggested an improved method for the scour depth prediction in clear water conditions, a well-designed physical model is recommended to investigate the scour characteristics under live bed conditions. In addition, different sediment sizes and non-uniform size sediment should be incorporated in the future research, as natural rivers generally consist of non-uniform sediment. Finally, scour modelling using computational fluid dynamics (CFD) should also be explored to refine the proposed equation for design purposes.

**Author Contributions:** S.H.H. and S.O.L. provided background data and motivation. S.H.H. and S.O.L. designed the experiments and determined the results; R.S. analyzed the data; R.S. and S.H.H. contributed to motivation, data analysis input, and interpretation of results; R.S. and S.H.H. interpreted the results and wrote the paper; and all authors participated in final review and editing of the paper.

**Funding:** National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (NRF-2017R1A2B2011990); West Virginia University internal grant.

**Acknowledgments:** The laboratory data used in this paper were measured in Georgia Tech and the permission was granted by Seung Ho Hong and Seung Oh Lee.

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

### **References**


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

### *Article* **Collapsing Mechanisms of the Typical Cohesive Riverbank along the Ningxia–Inner Mongolia Catchment**

### **Guosheng Duan 1, Anping Shu 1,\*, Matteo Rubinato 2, Shu Wang <sup>1</sup> and Fuyang Zhu <sup>1</sup>**


Received: 7 August 2018; Accepted: 11 September 2018; Published: 18 September 2018

**Abstract:** As one of the major sediment sources in rivers, bank collapse often occurs in the Ningxia–Inner Mongolia catchment and, to date, it caused substantial social, economic and environmental problems in both local areas and downstream locations. To provide a better understanding of this phenomenon, this study consisted of modifying the existing Bank Stability and Toe Erosion Model (BSTEM), commonly used to investigate similar phenomena, introducing new assumptions and demonstrating its applicability by comparing numerical results obtained against field data recorded at six gauging stations (Qingtongxia, Shizuishan, Bayan Gol, Sanhuhekou, Zhaojunfen, and Toudaoguai). Furthermore, the impact of multiple factors typical of flood and dry seasons on the collapse rate was investigated, and insights obtained should be taken into consideration when completing future projects of river adaptation and river restoration.

**Keywords:** Ningxia–Inner Mongolia; Yellow River; riverbank collapse; BSTEM model; flood & dry season; sediment transport

### **1. Introduction**

Riverbank collapse is a phenomenon caused by several natural factors (e.g., intense rainfall, site topography, properties of the river bed and hydraulic conditions of the river flow) [1] and artificial factors (e.g., man-made bank undercutting, basal clean-out) [2]. The Yellow River is the China's second-longest river and previous studies [3] confirmed that 518.38 km2 of the riverbank along the Ningxia–Inner Mongolia part of the Yellow River was eroded between 1958 and 2008. This elevated rate of riverbank erosion is due to continuous variations of flow levels, causing the deposition of loose riverbed material in alluvial plains, especially after heavier rainfall events typical of climate change [4,5]. To date, riverbank erosion has caused substantial social (e.g., traffic disruption due to flooding), economic (e.g., loss of farm land) and environmental problems (e.g., sediment dynamics and water quality) in both local areas and downstream locations [6,7].

Previous studies were completed to tackle this challenge and to identify what can induce riverbanks to collapse in catchments in both developing and developed countries. Grabowski [8] found that the soil texture, the soil structure, the unit weight and water composition [9,10] all play a principal role on the banks erosion. While non-cohesive banks are eroded through discrete particle entrainment that can be quantified using the magnitude of the shear stress and the particle size [11,12], cohesive banks are eroded through entrainment of aggregates [13] and this mechanism is very challenging to diagnose when considering the electrochemical forces acting between them [14,15]. The detachment and erosion phenomena of cohesive material by gravity and/or flowing water are typically controlled by a variety of physical, electrical and chemical forces, such as cohesion, electrochemical forces, pore-water pressure and matric suction [15]. These forces within and between aggregates cause their erosion to be complex. Furthermore, cohesive riverbanks are usually poorly drained due to their silt and clay composition, and thus may experience excess pore water pressures, one of the main agents of subaerial erosion. Clayey banks are also susceptible to desiccation cracking and slaking from wetting–drying cycles [16]. Numerous studies have examined the impact of water content on the riverbank erosion [17–22] while other studies have investigated the subaerial erosion due to the weakening and weathering of bank materials [23–25]. Additionally, Yu [26] completed a series of experimental tests to quantify the influence and the impact of the water infiltration on bank collapse and consequent sediment deposition. They identified that partial erosion can be due to fluvial hydraulic forces or gravitational forces. Moreover, total erosion of the bank can be generated if both forces are combined at their maximum magnitude, with a consequent disintegration of the bank and the transportation of fallen blocks within the flow, changing the bed elevation. Focusing on the Ningxia–Inner Mongolia banks, the erosion process has been found to be influenced by three specific factors: (i) the sand blown by the wind from the desert [27–30]; (ii) the sediment transported by the river from the upstream catchment [31,32]; and (iii) the material falling from the bank due to natural collapse [33]. It is fundamental to replicate all these governing processes within numerical models when assessing stability of the riverbanks. A large amount of literature is available and it covers all the numerical modeling aspects investigated to date [34–38]. One of the numerical tools most commonly used is the Bank Stability and Toe Erosion Model (BSTEM), developed by the National Sedimentation Laboratory in Oxford, Mississippi, USA [39], which has been continually modified and improved by the authors since its creation to improve its performance. To model the riverbanks' stability, BSTEM calculates a factor of safety (*Fs*) for multi-layers combining three equilibrium-method models: (i) the layers simulated are horizontal; (ii) vertical slices simulate tension crack; and (iii) the common failure mechanisms are cantilever failures. BSTEM also assumes that all the collapsed material is totally removed from the bank and does not accumulate on the toe [40]. BSTEM is frequently used to simulate: (i) the banks' stability and the consequent sediment loading inside the river [41]; (ii) stream rehabilitation projects [42]; and (iii) erosion and failure mechanisms [43,44]. Despite several findings that have resulted in a better understanding of this phenomenon, to the authors' knowledge there is a lack of studies which involve the evaluation of BSTEM for long-term streambank erosion and consequent long-term failure.

Therefore, to fill this gap, this work comprises a collection of field datasets in the Ningxia–Inner Mongolia area that were implemented on the BSTEM model to obtain a new solution to the mechanisms of long-term bank erosion. To achieve that, three hypothesis different from those applied on the original BSTEM model were executed: (i) when the soil collapses from the bank, part of it is washed away by the water, and the remaining part accumulates on the toe; (ii) the shape of the material that accumulates on the toe is considered to have a triangle form considering the cross-sectional view and the angle away from the river bank is the angle of repose of the sediment, which is taken as 30 degrees; and (iii) the downstream area of a bank is considered the first one to be affected by the collapse process.

The paper is organized as follows. In the next sections the study area, the procedure to collect field data and the numerical analysis applied to replicate the river bank collapse recorded (Section 2) are described. Section 3 presents the comparison between numerical and experimental datasets and the factors that could influence the rate of collapse are discussed. Finally, Section 4 reports the conclusions.

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

### *2.1. Case Study*

The Ningxia–Inner Mongolia catchment (Figure 1) is located in the lower part of the upper Yellow River and has a length of 913.5 km, starting from Qingtongxia (Ningxia Hui autonomous region) and ending at Hekouzhen (Inner Mongolia autonomous region). This site includes six gauging stations (Qingtongxia, Shizuishan, Bayan Gol, Sanhuhekou, Zhaojunfen, and Toudaoguai) and it is characterized by four desert regions (HedongSandyland, Ulan Buh Desert, Tengger Desert and Hobq Desert) and two alluvial plains (Yinchuan Plain and Hetao Plain).

The climate of this region can be classified as temperate continental climate, which is characterized by long cold winters, and warm to hot summers [45]. The annual rainfall (maximum and minimum peaks commonly reached in summer and winter, respectively) ranges from 100 to 288 mm.

The case study selected is located in Dengkou which is in the Inner Mongolia Province. The gauging station Bayan Gol is adjacent to it. Being part of the Hetao plain, the Dengkou catchment is characterized by a low gradient (≈0.17‰) and loose riverbed materials. This site was carefully chosen based on three specific criteria. Firstly, preliminary monitoring datasets of this site have confirmed that the riverbank was actively retreating. Secondly, the site was sufficiently close to a gauging station (about 5 km), hence detailed rainfall readings including time-series, frequencies and water surface elevations associated with each flow condition flowing next to the bank site were available and could be collected. Last but not least, the morphology and the soil properties of this area are representative of the cohesive streambank along all the Ningxia–Inner Mongolia catchment which has been drastically eroded from 1958 to 2008 [3]. Therefore, it is fundamental to provide a better understanding of this phenomenon, especially in this zone, to provide mitigation measures and plan actions to reduce the effects due to erosion of the riverbanks.

**Figure 1.** Location of the case study considered for this work. 1, 2, 3, 4 and 5 are the five surveying cross sections (500 m far from each other) selected for the collection of field-data.

### *2.2. Field Data*

To analyze the complex processes of riverbank collapse affecting Dengkou, it was necessary to inspect the nature of the bank and the typical flow conditions in the parts of the Yellow river crossing this area. Five cross sections (1, 2, 3, 4 and 5, Figure 1), spaced 500 m from each other, were identified along the banks for the acquisition of field datasets. Parameters and methods selected by the authors for the monitoring and collection of datasets from the field study are displayed in Table 1.


**Table 1.** Indexes and methods for in-site measurement.

Additionally, to measure the magnitude of bank collapse along the area of study, five stakes were positioned vertically to the flow direction at each section selected for the monitoring (Figure 2). Five measurements were conducted respectively on 21 May, 13 July, 9 August, 23 September and 30 September in 2011 to quantify the collapse distance. The collapse distances between each section recorded are displayed in Table 2.

**Table 2.** Collapse distances (m) recorded in 2011.


**Figure 2.** Collapse distance measurement for surveying Section 1.

### 2.2.1. Hydrologic Data

Monthly average hydrologic datasets (2011, year of site inspection) obtained from the gauging station Bayan Gol used for this study are showed in Table 3.


**Table 3.** Monthly average hydrologic data (2011).

### 2.2.2. Bank Characteristics of Each Section

Riverbank properties at each section (1, 2, 3, 4 and 5, Figure 1) were obtained through five field inspections in 2011 (1. 21 May, 2. 13 July, 3. 9 August, 4. 23 September and 5. 30 September), and details are showed in Table 4. During the characteristic long cold winters in this region, the Yellow River typically freezes, hence the summer season has been selected as sampling period to avoid the quantification of additional effects on river bed changes [46,47] which could have made the already complex comparison between experimental datasets and numerical results even more challenging.

**Table 4.** Characteristics of the sections selected for this study.


### *2.3. Riverbank Collapse Characterisation*

When collecting field data, observations were made to get more insights about the erosion that typically affects the riverbanks of the Ningxia–Inner Mongolia catchment. Figure 3 shows an example of the real case study under investigation prior and after the erosion. Observations confirmed that the bank material is cohesive, the riverbank slope is steep (77–86◦), and the pattern of bank collapse is considered to be planar, in conformity with previous studies conducted [48,49]. Additionally, while monitoring the site, it was possible to identify different stages of the continuous collapse process characterized as follows: (i) the bank toe erosion initiates (Figure 4a); (ii) tension cracks develop (Figure 4b); (iii) shearing starts and parts of the bank breaks and detaches from the main body (Figure 4c); (iv) bank failure occurs (Figure 4d). When the soil collapses from the bank, part of it is washed away by the water and the remaining part accumulates at the bank toe. The new shape of the bank including the additional material eroded then becomes the initial form of the next bank collapse.

**Figure 3.** Example of the existing bank slope (**left**) and an example of riverbank collapse along Dengkou reach (**right**).

**Figure 4.** Typical riverbank processes observed at the case study of Ningxia–Inner Mongolia catchment.

### *2.4. Numerical Modeling*

In this study, the BSTEM model, developed by the National Sedimentation Laboratory in Oxford, MS, USA, was used and modified according to the assumptions presented in the introduction.

### 2.4.1. BSTEM Method

There are two separate modules in the conventional BSTEM model: (i) the toe erosion module and (ii) the stability of the bank module [39]. Equation (1) is commonly applied to predict the width of the bank toe eroded due to the hydraulic conditions impacting on the cohesive riverbank [16,50–52]:

$$B = \kappa (\boldsymbol{\pi} - \boldsymbol{\pi}\_{\mathcal{E}})^{\boldsymbol{a}} \Delta t \tag{1}$$

where *B* is the toe erosion width (m); *κ* is the erodibility coefficient (m<sup>3</sup> N−<sup>1</sup> s−1); *τ* is the average shear stress (Pa); *τ<sup>c</sup>* is soil's critical shear stress (Pa); a is an exponent usually assumed to be unity, Δ*t* is the time interval (s). *κ* and *τ<sup>c</sup>* are functions of the soil properties and they are characterized by the following relationship, *κ* = 0.2*τ*−0.5 *<sup>c</sup>* . For non-cohesive soils, *τ<sup>c</sup>* is typically estimated based on the median particle diameter of the soil [53]. Rinaldi [54] noted the difficulty to accurately estimate *κ* but despite that, other methods provided solutions via a variety of methods to calculate *κ* and *τc*. One of

these methods was developed by Hanson [55] using an in situ jet-test device and his experimental results are currently used universally. In Equation (1), *τ* can be calculated using Equation (2):

$$
\pi = \gamma\_w RS \tag{2}
$$

where *γ<sup>w</sup>* is the weight of water (N/m3); *R* is the hydraulic radius (m); *S* is the channel slope.

For the specific scenario of planar failure assumed for this study, the riverbank stability analysis involved the computation of a safety factor, identified as the ratio between the resisting and the driving forces applied to the failure zone, which is typically calculated by using Equation (3) [56–59]:

$$F\_s = \frac{S\_R}{S\_D} \tag{3}$$

where *Fs* is the safety factor with respect to bank failure (*Fs* < 1 indicates the possibility of an eventual collapse); *SR* is the shear strength of the soil (kPa); *SD* is the driving stress (Pa). *SR* can be calculated using Equation (4):

$$S\_R = c' + S \tan \phi^b + \Psi \tan \phi' \tag{4}$$

where *c'* is the effective cohesion (the effective cohesion of the soil is the cohesion of the soil in the anti-shear process) (Pa); *S* is the normal stress (Pa); *φ<sup>b</sup>* is an angle that describes the relationship between the shear strength and the matric suction (◦); Ψ is the matric suction which contrasts the pore-water pressure (Pa); *φ*' is the effective internal friction angle (◦). *SD* can be calculated using Equation (5):

$$S\_D = \mathcal{W}\sin\beta\tag{5}$$

where *W* is the weight of the wet soil per unit area of the failure plane (Nm−2); *β* is the angle of the failure plane (◦) [39]. Note that the BSTEM model takes *c*' into account to represent the cohesion effect. Nevertheless, the influence of cohesive materials on erosion would be considered through other approaches, as well. For example, Dodaro et al. [60,61] modified Shields parameter to predict erosion phenomena at a cohesive sediment mixture.

### 2.4.2. BSTEM Method Modified

As defined in Section 2.3, the phenomenon of bank collapse is characterized by four different stages. Despite multiple factors influencing the rate of the riverbank collapse, as previously mentioned, continuous collapses are also affected by the intensity of the antecedent collapse magnitude. When the collapse happens, the riverbank affected typically divides into two sections within a moment-frame: one part collapsed is naturally washed away by the river and the remaining part usually deposits at the bank toe. Before the next bank collapse takes place, the additional deposited soil is eroded. Therefore, it is very important to characterize the continuous changing shape of the riverbank to correctly predict the effects of the continuous erosion. Although the original BSTEM model provides some insights into this process and is very helpful for the analysis of the riverbank stability, it should be noted that BSTEM assumes that the falling part which has collapsed from the riverbank does not accumulate on the toe and is immediately removed. Hence it does not accurately represent the reality of the entire erosion process typical on natural rivers. Authors observed multiple shapes and forms during the erosion process following the collapse of material and have tried to characterize each of them, but the literature published to date lacks formulae to accurately estimate and quantify the amount of material deposited. Therefore, to address this gap, three changes were made in the original BSTEM model as follows.

(i) The quantity of soil falling into the river and carried away by the natural streamflow in the form of suspended material can be calculated using Equation (6):

$$M\_0 = (S\_\* - S\_{\overline{v}})Hd\tag{6}$$

where *M*<sup>0</sup> is the weight of soil carried away per unit length *(*kg/m); *Sv* is the sediment concentration (kg/m3), typically obtained by field monitoring; *H* is the height of the bank (m); *d* is the width of the flow in the proximity of the bank (m). Because there was a main current which was about 20 m away from the bank, *d* is assumed as a constant (*d* = 20). *S*\* is the capacity of the flow to transport suspended material (kg/m3) and can be calculated using Equation (7) [62]:

$$\mathcal{S}\_{\*} = \mathbb{K} \left( \mathcal{U}\_{L}^{3} / \text{gh}\omega \right)^{\text{m}} \tag{7}$$

where *K* and *m* are parameters that can be obtained from the literature [60]; *UL* is flow velocity (m/s); *g* is gravitational acceleration (m/s2); *h* is water depth (m); and *ω* is the velocity of the particles setting (m/s) which can be calculated by using Equation (8) [62]:

$$
\omega = (\gamma\_{\text{s}} - \gamma\_{\text{w}}) \text{g} \text{D}^2 / 25.6 \gamma\_{\text{w}} \text{v} \tag{8}
$$

where *D* is the grain size of the material (m); *υ* is kinematic viscosity coefficient (m2/s).

(ii) The quantity of deposition soil can be calculated by using Equation (9):

$$M\_d = (W\_0 - M\_0)\mathbf{x} \tag{9}$$

where *Md* is the quantity of material deposited (kg); *W*<sup>0</sup> is the quantity of material collapsed (excluding the part suspended in the flow) (kg); *x* is a ratio between the material deposited and the material collapsed.

(iii) The new shape generated by the deposition of material on the bank toe is assumed to be triangular and the angle away from the riverbank is the angle of repose of sediment. In general, loose sediment grains in water accumulate with an angle close to the sediment angle of repose [63,64].

In the calculation process, the time step is selected as the hour. This decision is justified by the fact that since the flow rate, water depth, flow rate and other data are used as a daily average, the 24 h cycle was considered more appropriate. After the original boundary has passed 24 h, the river bank stability was tested: if the river bank was still stable, then the next set of water depth and flow rate conditions was inserted and the river bank toe erosion module continued to run. When the river bank stability coefficient was less than 1, the collapse occurred and consequently the collapse rate should have been recorded. For this scenario, the amount of collapse was calculated according to the shape of the collapse area identified, and the amount of sediment deposited at the foot of the slope was calculated by using Equations (6)–(9). The shape of the material deposited on the slope was considered to be triangular, and the angle away from the river bank was the angle of repose of the sediment, which was taken as 30 degrees. Based on the collapse rate, the new shape of the river bank was obtained and taken as a new initial condition for the next collapse.

The scheme utilized by the modified BSTEM model can be summarized with the flow chart displayed in Figure 5.

It is complex to quantify in real time the collapse processes along the entire riverbank of the study area due to flow conditions under continuous change, but it is possible to provide an approximation of the riverbank collapse processes by using averaged parameters such as monthly average water depths, flowrates, sediment concentration and average grain size. In addition, there are also parameters that require an accurate calibration during the entire calculation process, such as *x*. *τ<sup>c</sup>* can be obtained by using values furnished in the literature [65]. By applying trial and error techniques, *x* = 0.39 corresponds to the optimal numerical value to replicate the realistic distance of collapse measured for this case study and discrepancies obtained between numerical and field measurements are listed in Tables 5–9.

**Figure 5.** Flow chart of the established model.

**Table 5.** Numerical results, field measurements and errors for surveying Section 1 of Figure 1.


**Table 6.** Numerical results, field measurements and errors for surveying Section 2 of Figure 1.


**Table 7.** Numerical results, field measurements and errors for surveying Section 3 of Figure 1.



**Table 8.** Numerical results, field measurements and errors for surveying Section 4 of Figure 1.

**Table 9.** Numerical results, field measurements and errors for surveying Section 5 of Figure 1.


By analyzing the results displayed in Tables 5–9, it is possible to confirm that the error between numerical and experimental values is between 6% and 20%. There are multiple influencing factors that can affect this comparison and the shape of the bank toe, always changing under the action of water flow which provides the higher impact. As water and sand conditions are continuously varying, the error is considerably higher considering measurements vs numerical results for the smallest time frame monitored (23 September–30 September). On the other hand, for longer time scales, the total error is about 5%. Overall, the numerical results are fairly consistent with the data recorded during the monitoring campaign, considering the significant modifications applied within the numerical model.

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

This section provides the further results obtained and their interpretation.

### *3.1. Collapse Processes Obtained through Numerical Simulations with the Modified BSTEM Model*

The riverbank collapses simulated with the modified BSTEM model for the period 21 May–30 September in 2011 are illustrated in Figures 6–10 for surveying river Sections 1–5 respectively. Collapse rates of study sections in different months are listed in Table 10.

**Figure 6.** Collapse phases for surveying Section 1.

**Figure 10.** Collapse phases for Section 5.


**Table 10.** Collapse rates of study sections in different months (Unit: m/d).

To identify hypothetical relationships between hydraulic conditions and collapse rates, different indicators recorded were plotted versus the collapse rates estimated, and factors obtained are displayed in Figure 11.

**Figure 11.** The relationship between flow conditions and collapse rates, the effects of (**a**) flowrate; (**b**) water depth; (**c**) velocity; (**d**) sediment carrying capacity; and (**e**) bank height.

Because the flowrate, water depth, velocity and sediment carrying capacity were of similar magnitudes from June to July, the average collapse rate was used to discuss the effects on collapse rates by influencing factors previously highlighted, and their impacts can be observed in Figure 11.

The authors have attempted to characterize the continuous collapses influenced by multiple factors that can be divided into two categories: autologous and external. The autologous factor considered is the bank height, while external factors include flow rate, water depth, velocity and sediment carrying capacity. Different factors play dissimilar roles on influencing continuous collapse processes. For example, when flow rate, velocity and sediment carrying capacity rise, bank toe erosion can consequently increase due to higher shear stresses. Additionally, when water depth increases, soil shear strength diminishes due to the soil saturation, however the rise of water depth could also improve the bank stability due to the water pressure acting on the bank.

Even if bank collapses are continuously affected by these factors, there are still some interesting regular features observed in this study. For example, in Figure 11a–d, it can be observed that the collapse rate is positively correlated with the flow rate, water depth, velocity and capacity of carrying materials. Although the bank heights are all different in altered sections, collapse rates for the locations considered are higher in the flood season (July to September) than in the dry season (May to June), due to the increase of flow rates, water depth and velocities. Due to the strength of higher flow rates, the flow shear stress is typically greater than the soil shear stress, causing continuous erosion of the bank toe. It should also be noticed then that flow conditions are playing a leading role in the collapse processes and are the main dynamic factors for the quantity of material deposited.

Another interesting phenomenon which can be noticed within the results is that continuous collapse rates are bigger in May than in June, except for Section 3 (Figure 11). During this period, which belongs to the dry season, flow rates, water depth, velocities and sediment carrying capacities are all smaller (Table 3). This is due to the associated phenomena when water depths reduce, with a consequent loss in river bank stability and more intense bank collapses. However, for Section 3, the collapse distance decreases, and the authors believe that this is due to the bank height (2.3 m, the highest among the sections).

### *3.2. Non Applicability of the Traditional BSTEM Method*

In order to further demonstrate that the original BSTEM model is not applicable for the representation of continuous long term bank collapses, numerical and experimental results are listed in Tables 11–15, together with the errors obtained.


**Table 11.** BSTEM results, field measurements and errors for surveying Section 1.

**Table 12.** BSTEM results, field measurements and errors for surveying Section 2.



**Table 13.** BSTEM results, field measurements and errors for surveying Section 3.

**Table 14.** BSTEM results, field measurements and errors for surveying Section 4.


**Table 15.** BSTEM results, field measurements and errors for surveying Section 5.


It can be concluded that the errors are all very large for longer monitoring periods. That is due to the hypothesis implemented within the BSTEM method, commonly used to calculate collapse distance, considering no soil deposition at the bank toe after the collapse—an unrealistic condition.

### *3.3. Discussions*

Comparing simulation results using the modified BSTEM model with field data collected, as showed in Tables 5–9, it can be observed that errors range from 0.3% to 5.07%, hence are all acceptable and it is possible to confirm the reliability and accuracy of the modified BSTEM model to calculate continuous riverbank collapse.

Furthermore, additional outcomes have to be highlighted to increase the performance of the numerical model. Firstly, as the initial boundary conditions can have an impact on river bank collapse, supplementary studies have to be conducted to investigate this phenomenon. Secondly, due to the lack of advanced experimental measurement techniques, it is very difficult to get the river bed profile at every time step, and for this study a simplified bank shape was used in the calculations presented. This may be one more reason for the errors calculated. To solve this problem, sonar and electromagnetic wave measurement technologies may be adopted in future studies. Thirdly, in natural rivers, flow conditions such as flow rate and water depth vary continuously. To obtain bank collapse distances, monthly average flow conditions were used in this study, hence more targeted and accurate flow conditions should be used in future studies to improve the results. Finally, riverbed erosion and bank erosion typically occur at the same time and both are influenced by each other. For example, riverbed erosion influences the bank collapse by increasing the bank height and when the bank collapse happens, falling soil directly settles along the area where riverbed erosion happened before depositing at the bank toe.

As all these aspects are crucial, by implementing them within the numerical model it is predicted that results will be more accurate.

### **4. Conclusions**

(1) As confirmed by the observations made during the monitoring of the area, a single bank collapse can be divided into four different stages that can be summarized as: (i) bank toe erosion initiation; (ii) tension crack development; (iii) shearing initiation; and (iv) bank failure occurrence. When a single collapse happens, part of the collapsed soil is washed away by the water and the remaining part accumulates at the bank toe. The new shape of the bank formed after this process along the river becomes the initial form of the next bank collapse. Bank collapse processes in a longer time scale are continuous collapses that are made from single bank collapse during short time scales.

(2) Based on the original BSTEM model, a new modified version was implemented, applying several new assumptions. The original BSTEM assumed that the falling part which had collapsed from the riverbank did not accumulate on the toe and was immediately removed. This does not accurately represent the reality of the entire erosion process typical in natural rivers, and comparisons between experimental and numerical results have shown that, by using the original approach, significant errors will affect the numerical results, as listed in Tables 11–15.

Furthermore, to address this gap, a modified BSTEM was developed to obtain the quantity of collapsed sediment accumulating on the toe. Several parameters were introduced, such as sediment concentration, sediment carrying capacity and sediment setting velocity. The quantity of collapsed sediment accumulated on the toe could have been obtained using Formulas (6)–(9).

Additionally, another main goal modifying the BSTEM model was to obtain the new shape generated by the deposition of material on the bank toe. It should be noted that the boundary conditions play an important role in the calculations of continuous collapse processes for natural rivers. The new shape generated by the deposition of material on the bank toe was the initial boundary condition of next collapse process. Based on this, after observing multiple shapes during the erosion process following the collapse of material, the shape generated by the deposition of material on the bank toe was assumed to be triangular and the angle away from the riverbank was the angle of repose of sediment. After a single collapse, the new shape of the riverbank could have been obtained from the shape of the deposited sediment and the collapsed plain.

(3) Adopting all the previous modifications within the calculation of every single collapse process, the numerical results obtained using the modified BSTEM are consistent with the monitored datasets for continuous processes, as listed in Tables 5–9.

(4) The relationship between collapse rates and influencing factors was discussed and influencing factors were divided into two categories: autologous and external (flow conditions). Different factors play different roles in continuous collapse processes. For example, the flow rate, velocity and sediment carrying capacity rise can increase bank toe erosion by increasing flow shear stress. Water depth rise can reduce soil shear strength by promoting soil saturation, and water depth rise can also increase river bank stability by providing water pressure. Soil characteristics determine the shear strength of each material, and bank height plays an important role in the bank stability analysis. For similar flow conditions, a higher bank with a steeper slope is more unstable. These factors should be considered together to better characterize the collapse mechanisms.

(5) Aspects such as precisely determining the bank shape and quantifying the more singular flow conditions are crucial, and by implementing them within the numerical model it is predicted that results will become more accurate.

**Author Contributions:** All the authors jointly contributed to this research. A.S. was responsible for the proposition and design of the field monitoring and the modified BSTEM method; M.R. and G.D. analysed the experimental datasets and wrote the paper; and S.W. and F.Z. participated in the experiments.

**Funding:** The research reported in this manuscript is funded by the National Basic Research Program of China (Grant No. 2011CB403304) and National Natural Science Foundation of China (Grant No. 11372048).

**Acknowledgments:** The research reported in this manuscript is supported by the National Basic Research Program of China (Grant No. 2011CB403304) and National Natural Science Foundation of China (Grant No. 11372048).

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

### **References**


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